RE = function(n = 100000, hyper = c(0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0, 0, 100, 100), sessions = sessions){

n_m = length(unique(sessions[,2]))     #No. of months
n_s = dim(sessions)[2]-2 		#No. of sites (MSAs)
n_y = length(unique(sessions[,1]))	#No. of years

#Specify hyper-parameters using notation in calculations:
a = c(hyper[1], hyper[2])
b = c(hyper[3], hyper[4])
c = c(hyper[5], hyper[6])
d = c(hyper[7], hyper[8])
g = c(hyper[9], hyper[10])
h = c(hyper[11], hyper[12])

#Starting values for seasonal and site effects - site effects starting at their prior means, seasonal effects just starting at zero: 
start.gamma.mu = rep(0, n_m)
start.gamma.sigma = rep(0, n_m)
start.epsilon.mu = rep(0, n_s)
start.epsilon.sigma = rep(0, n_s)

#For trend:
x = seq(-1,1,length.out = n_y)
x.i = vector("numeric", length=dim(sessions)[1])
for(i in 1:length(x.i)){
	if(sessions[i,1]==2016){x.i[i]=x[1]}
	if(sessions[i,1]==2017){x.i[i]=x[2]}
	if(sessions[i,1]==2018){x.i[i]=x[3]}
	if(sessions[i,1]==2019){x.i[i]=x[4]}
	if(sessions[i,1]==2020){x.i[i]=x[5]}
	if(sessions[i,1]==2021){x.i[i]=x[6]}}

start.beta.mu = g[1]/h[1]
start.beta.sigma = g[2]/h[2]

#Storage for MCMC output:
X.mu = matrix(ncol = n_m, nrow = n)   #This matrix will take the MCMC output for the seasonal effects (gamma's) for mu - one row per iteration, one column per season
X.sigma = matrix(ncol = n_m, nrow = n)   #This matrix will take the MCMC output for the seasonal effects (gamma's) for (log)sigma - one row per iteration, one column per season

Y.mu = matrix(ncol = n_s, nrow = n)  #This matrix will take the MCMC output for the site effects (epsilon's) for mu - one row per iteration, one column per site
Y.sigma = matrix(ncol = n_s, nrow = n)  #This matrix will take the MCMC output for the site effects (epsilon's) for (log)sigma - one row per iteration, one column per site

Z.mu = vector("numeric",length = n)	#This vector will take the MCMC output for the trend effects for mu - one row per iteration
Z.sigma = vector("numeric",length = n)	#This vector will take the MCMC output for the trend effects for (log)sigma - one row per iteration


#Storage for candidates
can.X.mu = matrix(ncol = n_m, nrow = n)  #This matrix will take the candidates for the seasonal effects for mu
can.X.mu[1, ] = start.gamma.mu
can.X.sigma = matrix(ncol = n_m, nrow = n)  #This matrix will take the candidates for the seasonal effects for (log)sigma
can.X.sigma[1, ] = start.gamma.sigma

can.Y.mu = matrix(ncol = n_s, nrow = n)  #This matrix will take the candidates for the site effects for mu
can.Y.mu[1, ] = start.epsilon.mu
can.Y.sigma = matrix(ncol = n_s, nrow = n)  #This matrix will take the candidates for the site effects for (log)sigma
can.Y.sigma[1, ] = start.epsilon.sigma

can.Z.mu = vector("numeric", length = n)	#This vector will take the candidates for the trend effects for mu
can.Z.mu[1] = start.beta.mu
can.Z.sigma = vector("numeric", length = n)	#This vector will take the candidates for the trend effects for (log)sigma
can.Z.sigma[1] = start.beta.sigma

#Now we finalise the first row of the MCMC storage matrices at their initial states:
X.mu[1, ] = can.X.mu[1, ]
X.sigma[1, ] = can.X.sigma[1, ]
Y.mu[1, ] = can.Y.mu[1, ]
Y.sigma[1, ] = can.Y.sigma[1, ]
Z.mu[1] = can.Z.mu[1]
Z.sigma[1] = can.Z.sigma[1]

#Matrices for acceptance probabilities:
aprob.gamma.mu = matrix(ncol = n_m, nrow = n)
aprob.gamma.sigma = matrix(ncol = n_m, nrow = n)
aprob.epsilon.mu = matrix(ncol = n_s, nrow = n)
aprob.epsilon.sigma = matrix(ncol = n_s, nrow = n)
aprob.beta.mu = vector("numeric", length = n)
aprob.beta.sigma = vector("numeric", length = n)


#Log-likelihood function:
loglik = function(GAMMA.MU,GAMMA.SIGMA, EPSILON.MU,EPSILON.SIGMA, BETA.MU, BETA.SIGMA, DATA, X){
	N = length(DATA)
	bit1 = -(N/2)*log(2*pi)
	bit2 = sum(BETA.SIGMA*X+GAMMA.SIGMA+EPSILON.SIGMA)
	bit3 = sum(((DATA - BETA.MU*X - GAMMA.MU - EPSILON.MU)/(exp(BETA.SIGMA*X + GAMMA.SIGMA + EPSILON.SIGMA)))^2)
	llh = bit1 - bit2 -0.5*bit3
	return(llh)} 

#Initialise zeta.mu, zeta.sigma, tau.mu and tau.sigma using the starting values for the random effects:
zeta.mu = vector("numeric", n)
zeta.sigma = vector("numeric", n)
tau.mu = vector("numeric", n)
tau.sigma = vector("numeric", n)

zeta.mu[1] = rgamma(1, (n_s/2)+a[1], b[1]+0.5*sum(Y.mu[1, ]^2))
zeta.sigma[1] = rgamma(1, (n_s/2)+a[2], b[2]+0.5*sum(Y.sigma[1, ]^2))

grc = vector("numeric", n_m)
grc[1] = (X.mu[1,1]-0.5*(X.mu[1,n_m]+X.mu[1,2]))^2
grc[n_m] = (X.mu[1,n_m]-0.5*(X.mu[1,(n_m-1)]+X.mu[1,1]))^2
for(gg in 2:11){
	grc[gg] = (X.mu[1,gg]-0.5*(X.mu[1,(gg-1)]+X.mu[1,(gg+1)]))^2
	} 
tau.mu[1] = rgamma(1, (n_m/2)+c[1], d[1]+0.5*sum(grc))

grc2 = vector("numeric", n_m)
grc2[1] = (X.sigma[1,1]-0.5*(X.sigma[1,n_m]+X.sigma[1,2]))^2
grc2[n_m] = (X.sigma[1,n_m]-0.5*(X.sigma[1,(n_m-1)]+X.sigma[1,1]))^2
for(gg in 2:11){
	grc2[gg] = (X.sigma[1,gg]-0.5*(X.sigma[1,(gg-1)]+X.sigma[1,(gg+1)]))^2
	}
tau.sigma[1] = rgamma(1, (n_m/2)+c[2], d[2]+0.5*sum(grc2))

#Now start the MCMC/Gibbs iterations:

for(i in 2:n){  
	######################################### 
	########Seasonal effects for mu::########
        #########################################
	lik.X.mu = vector("numeric", n_m)
	lik.X.mu[] = 0
	for(m in 1:n_m){
		can.X.mu[i, m] = X.mu[(i-1), m] + rnorm(1, 0, err.gamma.mu[m])

	for(k in 1:n_s){
		data.contribution = sessions[, (k+2)][sessions[, 2]==m]
		x.contribution = x.i[sessions[,2]==m]

		lik.cont.X.mu = loglik(can.X.mu[i, m], X.sigma[(i-1),m], Y.mu[(i-1),k], Y.sigma[(i-1),k], Z.mu[i-1], Z.sigma[i-1], data.contribution, x.contribution) - loglik(X.mu[(i-1),m],X.sigma[(i-1),m],Y.mu[(i-1),k], Y.sigma[(i-1),k], Z.mu[i-1],Z.sigma[i-1], data.contribution, x.contribution)

		lik.X.mu[m] = lik.cont.X.mu + lik.X.mu[m]}
	m.before = m-1
	if(m.before==0){m.before=12}
	m.after = m+1
	if(m.after==13){m.after=1}

	aprob.gamma.mu[i, m] = min(1, (exp(lik.X.mu[m])*(dnorm(can.X.mu[i,m], mean(X.mu[(i-1), m.before],X.mu[(i-1), m.after]), 1/sqrt(tau.mu[i-1])))/(dnorm(X.mu[(i-1),m], mean(X.mu[(i-1), m.before], X.mu[(i-1), m.after]), 1/sqrt(tau.mu[i-1])))))

	u = runif(1)
	if(u < aprob.gamma.mu[i, m]){                                                 #################################
		X.mu[i, m] = can.X.mu[i, m]}else{X.mu[i, m] = X.mu[(i-1), m]}  }     #Seasonal effect for mu updated!#
        ############################################                                   #################################
	########Seasonal effects for sigma::########
        ############################################
	lik.X.sigma = vector("numeric", n_m)
	lik.X.sigma[] = 0
	for(m in 1:n_m){
		can.X.sigma[i, m] = X.sigma[(i-1), m] + rnorm(1, 0, err.gamma.sigma[m])

	for(k in 1:n_s){
		data.contribution = sessions[, (k+2)][sessions[, 2]==m]
		x.contribution = x.i[sessions[,2]==m]

		lik.cont.X.sigma = loglik(X.mu[i, m], can.X.sigma[i, m], Y.mu[(i-1),k], Y.sigma[(i-1),k], Z.mu[i-1], Z.sigma[i-1], data.contribution, x.contribution) - loglik(X.mu[i, m],X.sigma[(i-1),m],Y.mu[(i-1),k], Y.sigma[(i-1),k], Z.mu[i-1],Z.sigma[i-1], data.contribution, x.contribution)

		lik.X.sigma[m] = lik.cont.X.sigma + lik.X.sigma[m]}
	m.before = m-1
	if(m.before==0){m.before=12}
	m.after = m+1
	if(m.after==13){m.after=1}

	aprob.gamma.sigma[i, m] = min(1, (exp(lik.X.sigma[m])*(dnorm(can.X.sigma[i,m], mean(X.sigma[(i-1), m.before],X.sigma[(i-1), m.after]), 1/sqrt(tau.sigma[i-1])))/(dnorm(X.sigma[(i-1),m], mean(X.sigma[(i-1), m.before], X.sigma[(i-1), m.after]), 1/sqrt(tau.sigma[i-1])))))

	u = runif(1)
	if(u < aprob.gamma.sigma[i, m]){                                                          ########################################################
		X.sigma[i, m] = can.X.sigma[i, m]}else{X.sigma[i, m] = X.sigma[(i-1), m]}  }     #Now we have seasonal effects for mu and sigma updated!#
        #####################################                                                      ########################################################
	########Site effects for mu::########
        #####################################
	lik.Y.mu = vector("numeric", n_s) 
	lik.Y.mu[] = 0
	for(j in 1:n_s){
		can.Y.mu[i, j] = Y.mu[(i-1), j] + rnorm(1, 0, err.epsilon.mu[j])

	for(k in 1:n_m){
		data.contribution = sessions[, (j+2)][sessions[, 2]==k]
		x.contribution = x.i[sessions[,2]==k]

		lik.cont.Y.mu = loglik(X.mu[i, k], X.sigma[i, k], can.Y.mu[i, j], Y.sigma[(i-1),j], Z.mu[i-1], Z.sigma[i-1], data.contribution, x.contribution)-loglik(X.mu[i,k],X.sigma[i,k], Y.mu[(i-1),j], Y.sigma[(i-1),j], Z.mu[i-1], Z.sigma[i-1], data.contribution, x.contribution)
	
		lik.Y.mu[j] = lik.cont.Y.mu + lik.Y.mu[j]}

	aprob.epsilon.mu[i, j] = min(1, (exp(lik.Y.mu[j])*(dnorm(can.Y.mu[i,j], 0, 1/sqrt(zeta.mu[i-1])))/(dnorm(Y.mu[(i-1),j], 0, 1/sqrt(zeta.mu[i-1])))))

	u = runif(1)
	if(u < aprob.epsilon.mu[i, j]){                                            ##################################################################################################
		Y.mu[i, j] = can.Y.mu[i, j]}else{Y.mu[i, j] = Y.mu[(i-1), j]}}    #Now we have seasonal effects for mu and sigma updated, and also the site effects for mu updated!#
        ########################################                                    ##################################################################################################
	########Site effects for sigma::########
        ########################################
	lik.Y.sigma = vector("numeric", n_s)
	lik.Y.sigma[] = 0
	for(j in 1:n_s){
		can.Y.sigma[i, j] = Y.sigma[(i-1), j] + rnorm(1, 0, err.epsilon.sigma[j])

	for(k in 1:n_m){
		data.contribution = sessions[, (j+2)][sessions[, 2]==k]
		x.contribution = x.i[sessions[,2]==k]

		lik.cont.Y.sigma = loglik(X.mu[i, k], X.sigma[i, k], Y.mu[i, j], can.Y.sigma[i, j], Z.mu[i-1], Z.sigma[i-1], data.contribution, x.contribution)-loglik(X.mu[i,k], X.sigma[i,k], Y.mu[i, j], Y.sigma[(i-1),j], Z.mu[i-1], Z.sigma[i-1], data.contribution, x.contribution)
	
		lik.Y.sigma[j] = lik.cont.Y.sigma + lik.Y.sigma[j]}

	aprob.epsilon.sigma[i, j] = min(1, (exp(lik.Y.sigma[j])*(dnorm(can.Y.sigma[i,j], 0, 1/sqrt(zeta.sigma[i-1])))/(dnorm(Y.sigma[(i-1),j], 0, 1/sqrt(zeta.sigma[i-1])))))

	u = runif(1)
	if(u < aprob.epsilon.sigma[i, j]){                                                     ####################################################################################################
		Y.sigma[i, j] = can.Y.sigma[i, j]}else{Y.sigma[i, j] = Y.sigma[(i-1), j]}}    #Now we have seasonal effects for mu and sigma updated, and also the site effects for mu and sigma!#
        #####################################                                                   ####################################################################################################
	########Trend effect for mu::########
        #####################################
	can.Z.mu[i] = Z.mu[i-1] + rnorm(1, 0, err.beta.mu)  
	lik.Z.mu = 0
	for(j in 1:n_s){
		for(k in 1:n_m){
			data.contribution = sessions[, (j+2)][sessions[, 2]==k]
			x.contribution = x.i[sessions[,2]==k]
			lik.cont.Z.mu = loglik(X.mu[i, k], X.sigma[i, k], Y.mu[i, j], Y.sigma[i, j], can.Z.mu[i], Z.sigma[i-1], data.contribution, x.contribution) - loglik(X.mu[i,k], X.sigma[i,k], Y.mu[i,j], Y.sigma[i,j], Z.mu[i-1], Z.sigma[i-1], data.contribution, x.contribution)

	lik.Z.mu = lik.cont.Z.mu + lik.Z.mu}}

	aprob.beta.mu[i] = min(1, (exp(lik.Z.mu)*(dnorm(can.Z.mu[i], g[1], 1/sqrt(h[1])))/(dnorm(Z.mu[i-1], g[1], 1/sqrt(h[1])))))

	u = runif(1)
	if(u < aprob.beta.mu[i]){								###################################################################################################
		Z.mu[i] = can.Z.mu[i]}else{Z.mu[i] = Z.mu[i-1]}				#Now we have seasonal and site effects for mu and sigma updated. And also the trend effect for mu!#
        ########################################                                                ###################################################################################################
	########Trend effect for sigma::########
        ########################################
	can.Z.sigma[i] = Z.sigma[i-1] + rnorm(1, 0, err.beta.sigma)
	lik.Z.sigma = 0
	for(j in 1:n_s){
		for(k in 1:n_m){
			data.contribution = sessions[, (j+2)][sessions[, 2]==k]
			x.contribution = x.i[sessions[,2]==k]
			lik.cont.Z.sigma = loglik(X.mu[i, k], X.sigma[i, k], Y.mu[i, j], Y.sigma[i, j], Z.mu[i], can.Z.sigma[i], data.contribution, x.contribution) - loglik(X.mu[i,k], X.sigma[i,k], Y.mu[i,j], Y.sigma[i,j], Z.mu[i], Z.sigma[i-1], data.contribution, x.contribution)

	lik.Z.sigma = lik.cont.Z.sigma + lik.Z.sigma}}

	aprob.beta.sigma[i] = min(1, (exp(lik.Z.sigma)*(dnorm(can.Z.sigma[i], g[2], 1/sqrt(h[2])))/(dnorm(Z.sigma[i-1], g[2], 1/sqrt(h[2])))))

	u = runif(1)
	if(u < aprob.beta.sigma[i]){
		Z.sigma[i] = can.Z.sigma[i]}else{Z.sigma[i] = Z.sigma[i-1]} #ALL METROPOLIS UPDATES COMPLETE!													

	zeta.mu[i] = rgamma(1, (n_s/2)+a[1], b[1]+0.5*sum(Y.mu[i, ]^2))
	zeta.sigma[i] = rgamma(1, (n_s/2)+a[2], b[2]+0.5*sum(Y.sigma[i, ]^2))


##################################################################
grc.a = vector("numeric", n_m)
grc.a[1] = (X.mu[i,1]-0.5*(X.mu[i,n_m]+X.mu[i,2]))^2
grc.a[n_m] = (X.mu[i, n_m]-0.5*(X.mu[i, (n_m-1)]+X.mu[i,1]))^2
for(gg in 2:11){
	grc.a[gg] = (X.mu[i,gg]-0.5*(X.mu[i,(gg-1)]+X.mu[i,(gg+1)]))^2
	}
tau.mu[i] = rgamma(1, (n_m/2)+c[1], d[1]+0.5*sum(grc.a))

grc2.a = vector("numeric", n_m)
grc2.a[1] = (X.sigma[i,1]-0.5*(X.sigma[i,n_m]+X.sigma[i,2]))^2
grc2.a[n_m] = (X.sigma[i,n_m]-0.5*(X.sigma[i,(n_m-1)]+X.sigma[i,1]))^2
for(gg in 2:11){
	grc2.a[gg] = (X.sigma[i,gg]-0.5*(X.sigma[i,(gg-1)]+X.sigma[i,(gg+1)]))^2
	}
tau.sigma[i] = rgamma(1, (n_m/2)+c[2], d[2]+0.5*sum(grc2.a))			

###########################
#ALL UPDATES NOW COMPLETE!#
###########################
	print(i)
		}
	return(list(GAMMA.MU = X.mu, GAMMA.SIGMA = X.sigma, EPSILON.MU = Y.mu, EPSILON.SIGMA = Y.sigma, BETA.MU = Z.mu, BETA.SIGMA = Z.sigma, ZETA.MU = zeta.mu, ZETA.SIGMA = zeta.sigma, TAU.MU = tau.mu, TAU.SIGMA = tau.sigma, APROB.GAMMA.MU = aprob.gamma.mu, APROB.GAMMA.SIGMA = aprob.gamma.sigma, APROB.EPSILON.MU = aprob.epsilon.mu, APROB.EPSILON.SIGMA = aprob.epsilon.sigma, APROB.BETA.MU = aprob.beta.mu, APROB.BETA.SIGMA = aprob.beta.sigma))}