SSE = function(n = 100000, hyper = c(0.01,0.01,0.01, 0.01, 0,0.1,1000, 1000, 0.01, 0.01,0.01,0.01), data=months_data){
  
  n_m = 12 
  n_s = dim(data)[2]-1  
  
  a_sig = hyper[1]
  a_xi = hyper[2]
  b_sig = hyper[3]
  b_xi = hyper[4]
  d_sig = hyper[5]
  d_xi = hyper[6]
  e_sig = hyper[7]
  e_xi = hyper[8]
  f_sig = hyper[9]
  f_xi = hyper[10]
  g_sig = hyper[11]
  g_xi = hyper[12]
  
  start.gamma_sig = rep(0, n_m)
  start.epsilon_sig = rep(d_sig, n_s)
  start.gamma_xi = rep(0.1, n_m)
  start.epsilon_xi = rep(d_xi, n_s)
  
  err.gamma_sig = rep(0.05, n_m)
  err.gamma_xi = rep(0.1, n_m)
  err.epsilon_sig = rep(0.05, n_s)
  err.epsilon_xi = rep(0.1, n_s)
  
  W = matrix(ncol = n_m, nrow = n) #gamma_sigma
  X = matrix(ncol = n_m, nrow = n)   #This matrix will take the MCMC output for the seasonal effects (gamma's) for log(lambda) - one row per iteration, one column per season
  #gamma_xi
  Y = matrix(ncol = n_s, nrow = n)  #epsilon_sigma. This matrix will take the MCMC output for the site effects (epsilon's) for log(lambda) - one row per iteration, one column per site
  Z = matrix(ncol = n_s, nrow = n) #epsilon_xi
  
  can.W = matrix(ncol = n_m, nrow = n)  #This matrix will take the candidates for the seasonal effects
  can.W[1, ] = start.gamma_sig
  can.X = matrix(ncol = n_m, nrow = n)  #This matrix will take the candidates for the seasonal effects
  can.X[1, ] = start.gamma_xi
  
  can.Y = matrix(ncol = n_s, nrow = n)  #This matrix will take the candidates for the site effects
  can.Y[1, ] = start.epsilon_sig
  can.Z = matrix(ncol = n_s, nrow = n)  #This matrix will take the candidates for the site effects
  can.Z[1, ] = start.epsilon_xi
  
  #Now we finalise the first row of the MCMC matrices at their initial states:
  W[1, ] = can.W[1, ]
  X[1, ] = can.X[1, ]
  Y[1, ] = can.Y[1, ]
  Z[1, ] = can.Z[1, ]

  #Matrices for acceptance probabilities:
  aprob.gamma_sig = matrix(ncol = n_m, nrow = n)
  aprob.epsilon_sig = matrix(ncol = n_s, nrow = n)
  aprob.gamma_xi = matrix(ncol = n_m, nrow = n)
  aprob.epsilon_xi = matrix(ncol = n_s, nrow = n)
  
  #GPD log-likelihood function:
  loglik=function(GAMMA_sig,GAMMA_xi,EPSILON_sig,EPSILON_xi,dataset)
  {
    #ETA=GAMMA_sig+EPSILON_sig
    #XI=GAMMA_xi+EPSILON_xi
    m=min(1+(dataset*((GAMMA_xi+EPSILON_xi)/exp(GAMMA_sig+EPSILON_sig))))
    if(m<0.00001)return(as.double(-1000000))
    if(exp(GAMMA_sig+EPSILON_sig)<0.00001)return(as.double(-1000000))
    loglik=-length(dataset)*(GAMMA_sig+EPSILON_sig)-(1+(1/(GAMMA_xi+EPSILON_xi)))*sum(log(1+(((GAMMA_xi+EPSILON_xi)*dataset)/exp(GAMMA_sig+EPSILON_sig))))
    loglik
  }
  
  #Initialise tau, c and zeta using the starting values for the random effects:
  tau_sig = vector("numeric", n)
  c_sig = vector("numeric", n)
  zeta_sig = vector("numeric", n)
  tau_xi = vector("numeric", n)
  c_xi = vector("numeric", n)
  zeta_xi = vector("numeric", n)
  
  tau_sig[1] = rgamma(1, a_sig+(n_m/2), b_sig+0.5*sum(W[1, ]^2))
  c_sig[1] = rnorm(1, (d_sig*e_sig + (f_sig/g_sig)*sum(Y[1, ]))/(e_sig+n_s*(f_sig/g_sig)), 1/sqrt(e_sig+n_s*(f_sig/g_sig)))
  zeta_sig[1] = rgamma(1, f_sig+(n_s/2), g_sig+0.5*sum((Y[1, ]-c_sig[1])^2))
  tau_xi[1] = rgamma(1, a_xi+(n_m/2), b_xi+0.5*sum(X[1, ]^2))
  c_xi[1] = rnorm(1, (d_xi*e_xi + (f_xi/g_xi)*sum(Z[1, ]))/(e_xi+n_s*(f_xi/g_xi)), 1/sqrt(e_xi+n_s*(f_xi/g_xi)))
  zeta_xi[1] = rgamma(1, f_xi+(n_s/2), g_xi+0.5*sum((Z[1, ]-c_xi[1])^2))
  
  #Now start the MCMC/Gibbs iterations:
  
  for(i in 2:n){
    #Seasonal effects - sigma:
    lik.W = vector("numeric", n_m)
    lik.W[] = 0
    for(m in 1:n_m){
      can.W[i, m] = W[(i-1), m] + rnorm(1, 0, err.gamma_sig[m])
      
      for(k in 1:n_s){
        #Careful with you data - I'm using (k+1) here as the first column of my accident data was a month indicator (so no actual data)
        data.contribution = data[, (k+1)][data[, 1]==m]
        data.contribution = data.contribution[!is.na(data.contribution)]
        if(length(data.contribution)==0){lik.cont.W = 0}else{
          lik.cont.W = loglik(can.W[i, m], X[(i-1), m],Y[(i-1), k],Z[(i-1), k], data.contribution) - loglik(W[(i-1), m],  X[(i-1), m],Y[(i-1), k],Z[(i-1), k], data.contribution)}
        if(is.na(lik.cont.W)){print(k)}
        lik.W[m] = lik.cont.W + lik.W[m]}
      
      aprob.gamma_sig[i, m] = min(1, (exp(lik.W[m])*(dnorm(can.W[i, m], 0, 1/sqrt(tau_sig[i-1])))/(dnorm(W[(i-1), m], 0, 1/sqrt(tau_sig[i-1])))))
      
      u = runif(1)
      if(u < aprob.gamma_sig[i, m]){
        W[i, m] = can.W[i, m]}else{W[i, m] = W[(i-1), m]}  }  
    
    #Seasonal effects - xi:
    lik.X = vector("numeric", n_m)
    lik.X[] = 0
    for(m in 1:n_m){
      can.X[i, m] = X[(i-1), m] + rnorm(1, 0, err.gamma_xi[m])
      
      for(k in 1:n_s){
        #Careful with you data - I'm using (k+1) here as the first column of my accident data was a month indicator (so no actual data)
        data.contribution = data[, (k+1)][data[, 1]==m]
        data.contribution = data.contribution[!is.na(data.contribution)]
        if(length(data.contribution)==0){lik.cont.X = 0}else{
          lik.cont.X = loglik(W[i, m], can.X[i, m],Y[(i-1), k], Z[(i-1), k], data.contribution) - loglik(W[i, m], X[(i-1), m],Y[(i-1), k], Z[(i-1), k], data.contribution)}
        if(is.na(lik.cont.X)){print(k)}
        lik.X[m] = lik.cont.X + lik.X[m]}
      
      aprob.gamma_xi[i, m] = min(1, (exp(lik.X[m])*(dnorm(can.X[i, m], 0, 1/sqrt(tau_xi[i-1])))/(dnorm(X[(i-1), m], 0, 1/sqrt(tau_xi[i-1])))))
      
      u = runif(1)
      if(u < aprob.gamma_xi[i, m]){
        X[i, m] = can.X[i, m]}else{X[i, m] = X[(i-1), m]}  }  
    
    #Site effects-sigma:
    lik.Y = vector("numeric", n_s)
    lik.Y[] = 0
    for(j in 1:n_s){
      can.Y[i, j] = Y[(i-1), j] + rnorm(1, 0, err.epsilon_sig[j])
      
      for(k in 1:n_m){
        data.contribution = data[, (j+1)][data[, 1]==k]
        data.contribution = data.contribution[!is.na(data.contribution)]
        if(length(data.contribution)==0){lik.cont.Y = 0}else{
          lik.cont.Y = loglik(W[i, k], X[i, k],can.Y[i, j], Z[(i-1),j], data.contribution) - loglik(W[i, k], X[i, k],Y[(i-1), j], Z[(i-1),j], data.contribution)}
        lik.Y[j] = lik.cont.Y + lik.Y[j]}
      
      aprob.epsilon_sig[i, j] = min(1, (exp(lik.Y[j])*(dnorm(can.Y[i, j], c_sig[i-1], 1/sqrt(zeta_sig[i-1])))/(dnorm(Y[(i-1), j], c_sig[i-1], 1/sqrt(zeta_sig[i-1])))))
      
      u = runif(1)
      if(u < aprob.epsilon_sig[i, j]){
        Y[i, j] = can.Y[i, j]}else{Y[i, j] = Y[(i-1), j]}   }
    
    #Site effects-xi:
    lik.Z = vector("numeric", n_s)
    lik.Z[] = 0
    for(j in 1:n_s){
      can.Z[i, j] = Z[(i-1), j] + rnorm(1, 0, err.epsilon_xi[j])
      
      for(k in 1:n_m){
        data.contribution = data[, (j+1)][data[, 1]==k]
        data.contribution = data.contribution[!is.na(data.contribution)]
        if(length(data.contribution)==0){lik.cont.Z = 0}else{
          lik.cont.Z = loglik(W[i, k], X[i, k],Y[i, j], can.Z[i, j], data.contribution) - loglik(W[i, k], X[i, k],Y[i, j], Z[(i-1), j], data.contribution)}
        lik.Z[j] = lik.cont.Z + lik.Z[j]}
      
      aprob.epsilon_xi[i, j] = min(1, (exp(lik.Z[j])*(dnorm(can.Z[i, j], c_xi[i-1], 1/sqrt(zeta_xi[i-1])))/(dnorm(Z[(i-1), j], c_xi[i-1], 1/sqrt(zeta_xi[i-1])))))
      
      u = runif(1)
      if(u < aprob.epsilon_xi[i, j]){
        Z[i, j] = can.Z[i, j]}else{Z[i, j] = Z[(i-1), j]}   }
    
    tau_sig[i] = rgamma(1, a_sig+(n_m/2), b_sig+0.5*sum(W[i, ]^2))
    c_sig[i] = rnorm(1, (d_sig*e_sig + zeta_sig[i-1]*sum(Y[i, ]))/(e_sig+n_s*zeta_sig[i-1]), 1/sqrt(e_sig+n_s*zeta_sig[i-1]))
    zeta_sig[i] = rgamma(1, f_sig+(n_s/2), g_sig+0.5*sum((Y[i, ]-c_sig[i])^2))  #### All updates now complete!  
    tau_xi[i] = rgamma(1, a_xi+(n_m/2), b_xi+0.5*sum(X[i, ]^2))
    c_xi[i] = rnorm(1, (d_xi*e_xi + zeta_xi[i-1]*sum(Z[i, ]))/(e_xi+n_s*zeta_xi[i-1]), 1/sqrt(e_xi+n_s*zeta_xi[i-1]))
    zeta_xi[i] = rgamma(1, f_xi+(n_s/2), g_xi+0.5*sum((Z[i, ]-c_xi[i])^2))  #### All updates now complete! 
    
    print(i)
  }
  
  return(list(GAMMA_sig = W,GAMMA_xi = X, EPSILON_sig = Y,EPSILON_xi = Z, TAU_sig = tau_sig,TAU_xi = tau_xi, C_sig = c_sig,C_xi = c_xi, ZETA_sig = zeta_sig,ZETA_xi = zeta_xi, APROB.GAMMA_sig = aprob.gamma_sig,APROB.GAMMA_xi = aprob.gamma_xi, APROB.EPSILON_sig = aprob.epsilon_sig,APROB.EPSILON_xi = aprob.epsilon_xi))  
}


###Now run the code:
months_data = read.table("data_above_thr_all.txt")

test2=SSE(n = 50000, data=months_data)

print(c(mean(test$APROB.GAMMA_sig, na.rm=TRUE), mean(test$APROB.GAMMA_xi, na.rm=TRUE), mean(test$APROB.EPSILON_sig, na.rm=TRUE), mean(test$APROB.EPSILON_xi, na.rm=TRUE)))

par(mfrow=c(3,4))
for(i in 1:12){
  plot(ts(test2$GAMMA_sig[,i]))
}
for(i in 1:12){
  plot(ts(test2$GAMMA_xi[,i]))
}
for(i in 1:12){
  plot(ts(test2$EPSILON_sig[,i]))
}
for(i in 1:12){
  plot(ts(test2$EPSILON_xi[,i]))
}