bayes6<-function(it, dataset, ti, U, beta0start,beta1start,xistart,alphastart, errbeta0,errbeta1,errxi,erralpha, sdbeta0,sdbeta1,sdxi, lik.tol){ ###Notes on function arguments:: #1. "it" is the number of MCMC iterations #2. The input dataset - "dataset" - needs to be a two column vector for margins X and Y. In a time series context for temporal dependence modelling we would have X=X_t and Y=X_t+1 #3. "ti" is the time indicator (0/1 for before/after in the Canada application). The values in "ti" need to matched to the pairs (rows) in "dataset" so obviously there'll be one pair for which "ti" is both 0 and 1 - remove that pair? #4. "U" is the 2-D threshold for (X,Y). In a time series context, U[1]=U[2] #5. "beta0start" and "beta1start" each take two values, as does "xistart" - we have two margins (X and Y), which is not so obvious when fitting this model in a time series context! #6. Similarly, "errbeta0", "errbeta1" and "errxi" are vectors of length 2 - the random walk variances in the MCMC #7. And of course, "sdbeta0", "sdbeta1" and "sdxi" are also vectors of length 2 - the prior standard deviations in our independent N(0, sd*) priors #8. "alphastart" and "erralpha" are scalars, as we have a single parameter governing the extremal dependence (the logistic dependence parameter, alpha); we don't specify the prior standard deviation for alpha (would have been "sdalpha") as we are using a U(0,1) prior #Length of the dataset (NB: dataset is a 2-column matrix) - both X and Y need to be the same length, regardless of the application. Not joint threshold excesses, full pairwise observartions:: k = length(dataset[,1]) #Not actually sure if this is used anywhere... ###Vectors to take candidate values:: #First the marginal intercept and slope for log(sigma)_X - beta0(X) and beta1(X):: caneta_x0 = vector("numeric"); caneta_x0[1] = beta0start[1] #Intercept for X-margin caneta_x1 = vector("numeric"); caneta_x1[1] = beta1start[1] #Slope for X-margin #Now, the marginal intercept and slope for log(sigma)_Y - beta0(Y) and beta1(Y):: caneta_y0 = vector("numeric"); caneta_y0[1] = beta0start[2] #Intercept for Y-margin caneta_y1 = vector("numeric"); caneta_y1[1] = beta1start[2] #Slope for Y-margin #Next: the marginal shape parameters xi_X and xi_Y:: canxi_x = vector("numeric"); canxi_x[1] = xistart[1] canxi_y = vector("numeric"); canxi_y[1] = xistart[2] #Finally, the logistic dependence parameter alpha: canalpha = vector("numeric"); canalpha[1] = alphastart ###Vectors to take final (accepted) values in the MCMC scheme:: x0_eta = vector("numeric"); x0_eta[1] = caneta_x0[1] #Intercept for log(sigma)_X x1_eta = vector("numeric"); x1_eta[1] = caneta_x1[1] #Slope for log(sigma)_X y0_eta = vector("numeric"); y0_eta[1] = caneta_y0[1] #Intercept for log(sigma)_Y y1_eta = vector("numeric"); y1_eta[1] = caneta_y1[1] #Slope for log(sigma)_Y x_xi = vector("numeric"); x_xi[1] = canxi_x[1] #Shape xi_X y_xi = vector("numeric"); y_xi[1] = canxi_y[1] #Shape xi_Y z = vector("numeric"); z[1] = canalpha[1] #Logistic dependence parameter alpha ###Vectors to take the acceptance probabilities:: aprobeta_x0 = vector("numeric") #Intercept for log(sigma)_X aprobeta_x1 = vector("numeric") #Slope for log(sigma)_X aprobeta_y0 = vector("numeric") #Intercept for log(sigma)_Y aprobeta_y1 = vector("numeric") #Slope for log(sigma)_Y aprobxi_x = vector("numeric") #Shape xi_X aprobxi_y = vector("numeric") #Shape xi_Y aprobalpha = vector("numeric") #Logistic dependence parameter alpha ###Now the function which returns the log-likelihood for the bivariate threshold excess model (BTEM) with logistic dependence:: loglik = function(data, thresh, TI, ETA_x0,ETA_x1, ETA_y0,ETA_y1, XI_x, XI_y, ALPHA, LIK.TOL) ###Notes on arguments:: #1. "data" is, once again, a 2-column vector for margins X and Y (or X_t and X_t+1 in a time series context). We use "data" and not "dataset" here to keep the "loglik" function portable. data=dataset will be passed through on execution in the block updates below (see for loop) #2. "thresh" is, once again, a 2-D threshold (and of course thresh[1]=thresh[2] in a time series context). Again, a diferent label is used here to keep this function portable, although thresh=U will be passed through on execution in the block updates below #3. "TI" is the time indicator for each pair (row) in "data" (0/1 respectively for the before/after periods). TI=ti will be passed through on execution below #4. ETA_x0 is the value of the intercept for log(sigma)_X at which the BTEM likelihood is evaluated #5. ETA_x1 is the value of the slope for log(sigma)_X at which the BTEM likelihood is evaluated #6. ETA_y0 is the value of the intercept for log(sigma)_Y at which the BTEM likelihood is evaluated #7. ETA_y1 is the value of the slope for log(sigma)_Y at which the BTEM likelihood is evaluated #8. Xi_x is the value of the shape xi_X at which the BTEM likelihood is evaluated #9. Xi_x is the value of the shape xi_Y at which the BTEM likelihood is evaluated #10. ALPHA is the value of the logistic dependence parameter at which the BTEM likelihood is evaluated { u_x = thresh[1] #Marginal threshold label for X u_y = thresh[2] #Marginal threshold label for Y if (ALPHA<0.015 | ALPHA>0.99)return(-1000000) #Avoid numerical instabilities for very small or very large values of alpha (rarely seen in practice anyway) n = dim(data)[1] #Number of pairs lambda_x = sum(data[,1] > u_x)/n #Threshold excess rate for X margin lambda_y = sum(data[,2] > u_y)/n #Threshold excess rate for Y margin quad1 = data[data[,1] > u_x & data[,2] > u_y,] #Quadrant 1: Both observations in pair exceed their marginal thresholds quad2 = data[data[,1] > u_x & data[,2] <= u_y,] #Quadrant 2: Censoring: Only the observation in the X margin exceeds its threshold quad3 = data[data[,1] <= u_x & data[,2] > u_y,] #Quadrant 3: Censoring: Only the observation in the Y margin exceeds its threshold quad4 = data[data[,1] <= u_x & data[,2] <= u_y,] #Quadrant 4: Censoring: Both observations in pair lie sub-threshold #Times in quadrants. To model the dependence of log(sigma) on time, we need to corresponding timestamps (and hence values of "TI") for all pairs in each quadrant:: indquad1 = vector("numeric",length = dim(quad1)[1]) #Will take timestamps for each pair in quadrant 1 indquad2 = vector("numeric",length = dim(quad2)[1]) #Will take timestamps for each pair in quadrant 2 indquad3 = vector("numeric",length = dim(quad3)[1]) #Will take timestamps for each pair in quadrant 3 indquad4 = vector("numeric",length = dim(quad4)[1]) #Will take timestamps for each pair in quadrant 4 for(ii in 1:dim(data)[1]){ for(jj in 1:dim(quad1)[1]){ if(data[ii,]==quad1[jj,]){indquad1[jj]=ii}} #Record the timestamps for pairs in quadrant 1 (in the vector "indquad1") for(kk in 1:dim(quad2)[1]){ if(data[ii,]==quad2[kk,]){indquad2[kk]=ii}} #Record the timestamps for pairs in quadrant 2 (in the vector "indquad2") for(ll in 1:dim(quad3)[1]){ if(data[ii,]==quad3[ll,]){indquad3[ll]=ii}} #Record the timestamps for pairs in quadrant 3 (in the vector "indquad3") for(mm in 1:dim(quad4)[1]){ if(data[ii,]==quad4[mm,]){indquad4[mm]=ii}}} #Record the timestamps for pairs in quadrant 4 (in the vector "indquad4") tiquad1 = TI[indquad1] #Now record the corresponding values of "TI" for each pair in quadrant 1 tiquad2 = TI[indquad2] #Now record the corresponding values of "TI" for each pair in quadrant 2 tiquad3 = TI[indquad3] #Now record the corresponding values of "TI" for each pair in quadrant 3 tiquad4 = TI[indquad4] #Now record the corresponding values of "TI" for each pair in quadrant 4 ###Now work on evaluation of the BTEM likelihood; see (for example) Fawcett and Green (2018). Notice that, throughout, sigma = exp(beta0+beta1*t) to capture the treatment effect in the Canada analysis:: core1x = 1 + XI_x*(quad1[,1]-u_x)/exp(ETA_x0+(sum(tiquad1))*ETA_x1) core1y = 1 + XI_y*(quad1[,2]-u_y)/exp(ETA_y0+(sum(tiquad1))*ETA_y1) core2x = 1 + XI_x*(quad2[,1]-u_x)/exp(ETA_x0+(sum(tiquad2))*ETA_x1) core3y = 1 + XI_y*(quad3[,2]-u_y)/exp(ETA_y0+(sum(tiquad3))*ETA_y1) if(min(core1x)<0.00001)return(-1000000) if(min(core1y)<0.00001)return(-1000000) if(min(core2x)<0.00001)return(-1000000) if(min(core3y)<0.00001)return(-1000000) xtil1.a = (1+XI_x*(quad1[,1]-u_x)/exp( ETA_x0+(sum(tiquad1))*ETA_x1 ))**(-1/XI_x) xtil1.a[which(xtil1.a < 10**(-LIK.TOL-1))] = 10**(-LIK.TOL-1) xtil1=-(log(1-lambda_x * xtil1.a))**(-1) ytil1.a = (1+XI_y*(quad1[,2]-u_y)/exp( ETA_y0+(sum(tiquad1))*ETA_y1 ))**(-1/XI_y) ytil1.a[which(ytil1.a < 10**(-LIK.TOL-1))] = 10**(-LIK.TOL-1) ytil1=-(log(1-lambda_y * ytil1.a))**(-1) dxtil1bydx = ((lambda_x/exp( ETA_x0+(sum(tiquad1))*ETA_x1 ))*core1x**((-1/XI_x)-1))*((log(1-lambda_x*core1x**(-1/XI_x)))**(-2))/(1-lambda_x*core1x**(-1/XI_x)) dxtil1bydx[which(dxtil1bydx==Inf)]=1e+6 dytil1bydy = ((lambda_y/exp( ETA_y0+(sum(tiquad1))*ETA_y1 ))*core1y**((-1/XI_y)-1))*((log(1-lambda_y*core1y**(-1/XI_y)))**(-2))/(1-lambda_y*core1y**(-1/XI_y)) dytil1bydy[which(dytil1bydy==Inf)]=1e+6 bit1 = xtil1**(-1/ALPHA) + ytil1**(-1/ALPHA) loglik1 = -bit1**ALPHA + ((-1/ALPHA)-1)*(log(xtil1) + log(ytil1)) + (ALPHA-2)*log(bit1) +log(bit1**ALPHA -1 +1/ALPHA) + log(dxtil1bydx) + log(dytil1bydy) xtil2.a = (1+XI_x*(quad2[,1]-u_x)/exp( ETA_x0+(sum(tiquad2))*ETA_x1 ))**(-1/XI_x) xtil2.a[which(xtil2.a < 10**(-LIK.TOL-1))] = 10**(-LIK.TOL-1) xtil2=-(log(1-lambda_x * xtil2.a))**(-1) ytil2 = -log(1-lambda_y)**(-1) dxtil2bydx = ((lambda_x/exp( ETA_x0+(sum(tiquad2))*ETA_x1 ))*core2x**((-1/XI_x)-1))*((log(1-lambda_x*core2x**(-1/XI_x)))**(-2))/(1-lambda_x*core2x**(-1/XI_x)) dxtil2bydx[which(dxtil2bydx==Inf)]=1e+6 bit2 = xtil2**(-1/ALPHA) + ytil2**(-1/ALPHA) loglik2 = -bit2**ALPHA + ((-1/ALPHA) - 1)*log(xtil2) + (ALPHA-1)*log(bit2) + log(dxtil2bydx) xtil3 = -log(1-lambda_x)**(-1) ytil3.a = (1+XI_y*(quad3[,2]-u_y)/exp( ETA_y0+(sum(tiquad3))*ETA_y1 ))**(-1/XI_y) ytil3.a[which(ytil3.a < 10**(-LIK.TOL-1))] = 10**(-LIK.TOL-1) ytil3=-(log(1-lambda_y * ytil3.a))**(-1) dytil3bydy = ((lambda_y/exp( ETA_y0+(sum(tiquad3))*ETA_y1 ))*core3y**((-1/XI_y)-1))*((log(1-lambda_y*core3y**(-1/XI_y)))**(-2))/(1-lambda_y*core3y**(-1/XI_y)) dytil3bydy[which(dytil3bydy==Inf)]=1e+6 bit3 = xtil3**(-1/ALPHA) + ytil3**(-1/ALPHA) loglik3 = - bit3**ALPHA + ((-1/ALPHA) - 1)*log(ytil3) + (ALPHA-1)*log(bit3) + log(dytil3bydy) xtil4 = -log(1-lambda_x)**(-1) ytil4 = -log(1-lambda_y)**(-1) loglik4 = length(quad4[,1])*exp(-(((xtil4^(-1/ALPHA))+(ytil4^(-1/ALPHA)))^(ALPHA))) #llm = (sum(loglik1)+sum(loglik2)+sum(loglik3)+loglik4) #Naughty adjustment for now, to avoid the (very) occasional NaN's in the log-likelihood contributions for certain combinations of parameter values. #Where we get a NaN as the log-likelihood contribution for an individual pair, we replace it with the mean of the other log-likelihood contributions #Need to sort, but seems to do the job to keep it ticking over, and it only seems to happen very occasionally:: if(is.na(sum(loglik1))){S.loglik1=mean(loglik1, na.rm=TRUE)*length(loglik1)}else{S.loglik1=sum(loglik1)} if(is.na(sum(loglik2))){S.loglik2=mean(loglik2, na.rm=TRUE)*length(loglik2)}else{S.loglik2=sum(loglik2)} if(is.na(sum(loglik3))){S.loglik3=mean(loglik3, na.rm=TRUE)*length(loglik3)}else{S.loglik3=sum(loglik3)} llm = (S.loglik1+S.loglik2+S.loglik3+loglik4) return(llm) } for(i in 2:it){ print(i) #Print the MCMC iteration counter ######################FIRST BLOCK: UPDATE beta0_x ###################### caneta_x0[i] = x0_eta[i-1]+rnorm(1,0,errbeta0[1]) l1 = loglik(dataset, U, ti, caneta_x0[i], x1_eta[i-1], y0_eta[i-1], y1_eta[i-1], x_xi[i-1], y_xi[i-1], z[i-1], lik.tol)-loglik(dataset, U, ti, x0_eta[i-1], x1_eta[i-1], y0_eta[i-1], y1_eta[i-1], x_xi[i-1], y_xi[i-1], z[i-1], lik.tol) likely1 = exp(l1) aprobeta_x0[i] = min(1, likely1*dnorm(caneta_x0[i], 0, sdbeta0[1])/dnorm(x0_eta[i-1], 0, sdbeta0[1])) print(c(caneta_x0[i],x1_eta[i-1], y0_eta[i-1], y1_eta[i-1], x_xi[i-1], y_xi[i-1], z[i-1], aprobeta_x0[i])) #print(c("beta0(X)",aprobeta_x0[i])) u = runif(1) if(ufx) { x2<-0 while(x2