#x is a three dimensional array. The first argument of x denotes the panels; the second argument denotes time, the third argument denotes the functional component. #standard deviation function across panels std<-function(x){ return(apply(x,c(1,3),sd)) } #multiplier bootstrap boots<-function(x,m=floor((dim(x)[2])^(1/3)),stdx,repli=1000){ #r=dim(x)[1] n=dim(x)[2] #fres=dim(x)[3] cum=aperm(apply(x,c(1,3),cumsum),c(2,1,3)) psum=(cum[,(m+1):n,]-cum[,1:(n-m),])/m aa=cum[,n,]/n rm(cum) tsum=psum for(i in 1:(n-m)) tsum[,i,]=aa distr=rep(0,repli) psum=psum-tsum rm(tsum) #stdx=std(x) for(i in 1:repli){ g=rnorm(n-m) boo=aperm(aperm(psum,c(2,1,3))*g,c(2,1,3)) inter=apply(boo,c(1,3),sum)*sqrt(m/(n-m)) distr[i]=max(abs(inter)/stdx) } rm(psum) return(distr) } #joint simultaneous confidence bands jscb<-function(x,m=floor((dim(x)[2])^(1/3)),repli=1000,alpha=0.05){ rr=dim(x)[1] n=dim(x)[2] fres=dim(x)[3] k=length(alpha) stdx1=std(x) r=sort(boots(x,m,stdx=stdx1,repli)) crit=r[floor(repli*(1-alpha))] #upper=matrix(0,nrow=rr,ncol=fres) #lower=matrix(0,nrow=rr,ncol=fres) upper=array(0,dim=c(rr,fres,k)) lower=array(0,dim=c(rr,fres,k)) mm=apply(x,c(1,3),mean) for(i in 1:k){ upper[,,i]=mm+stdx1*crit[i]/sqrt(n) lower[,,i]=mm-stdx1*crit[i]/sqrt(n) } rm(stdx1) list(upper=upper,lower=lower,fmean=mm,crit=crit) } #test of parallelism Fptest<-function(x,m=floor((dim(x)[2])^(1/3)),repli=1000,alpha=0.05){ require(abind) r=dim(x)[1] n=dim(x)[2] fres=dim(x)[3] cx<-apply(x,c(1,2),mean) transx=x for(i in 1:fres) transx[,,i]=x[,,i]-cx rm(cx) xx<-array(transx[r,,]-transx[r-1,,],dim=c(1,n,fres)) for(i in (r-2):1){ for(j in r:(i+1)) xx<-abind(xx,transx[j,,]-transx[i,,],along=1) } rm(transx) stdx1=std(xx) r=sort(boots(xx,m,stdx=stdx1,repli)) mm=apply(xx,c(1,3),mean) rm(xx) crit=r[floor(repli*(1-alpha))] Toverall=sqrt(n)*max(abs(mm)/stdx1) Tpair=sqrt(n)*apply(abs(mm)/stdx1,1,max) rm(mm) rm(stdx1) list(Toverall=Toverall,Tpair=Tpair,crit=crit) } #minimum volatility method mv<-function(x,mrange=floor(dim(x)[2]^(1/3)/2):2*floor(dim(x)[2]^(1/3))){ r=dim(x)[1] n=dim(x)[2] fres=dim(x)[3] ddlrv=array(0,dim=c(r,fres,length(mrange))) cum=aperm(apply(x,c(1,3),cumsum),c(2,1,3)) for(i in 1:length(mrange)){ m=mrange[i] psum=(cum[,(m+1):n,]-cum[,1:(n-m),])/m aa=cum[,n,]/n tsum=psum for(j in 1:(n-m)) tsum[,j,]=aa psum=psum-tsum rm(tsum) ddlrv[,,i]=apply(psum^2,c(1,3),sum)*m/(n-m) } rm(cum) gg=rep(0,length(mrange)-2) for(i in 2:(length(mrange)-1)){ ggsd<-apply(ddlrv[,,(i-1):(i+1)],c(1,2),sd) gg[i-1]=mean(ggsd) } rm(ddlrv) list(vseries=gg,m=mrange[which(gg==min(gg))+1]) } ################################### ##generate Legendre polynomials ( for simulations) leg1<-function(x) return(rep(1,length(x))) leg2<-function(x) return(sqrt(3)*(2*x-1)) leg3<-function(x){ y=2*x-1 return(sqrt(5)*(3*y^2-1)/2) } leg4<-function(x){ y=2*x-1 return(sqrt(7)*(5*y^3-3*y)/2) } leg5<-function(x){ y=2*x-1 return(sqrt(9)*(35*y^4-30*y^2+3)/8) } leg6<-function(x){ y=2*x-1 return(sqrt(11)*(63*y^5-70*y^3+15*y)/8) } leg7<-function(x){ y=2*x-1 return(sqrt(13)*(231*y^6-315*y^4+105*y^2-5)/16) } leg8<-function(x){ y=2*x-1 return(sqrt(15)*(429*y^7-693*y^5+315*y^3-35*y)/16) } leg9<-function(x){ y=2*x-1 return(sqrt(17)*(6435*y^8-12012*y^6+6930*y^4-1260*y^2+35)/128) } leg10<-function(x){ y=2*x-1 return(sqrt(19)*(12155*y^9-25740*y^7+18018*y^5-4620*y^3+315*y)/128) } leg<-function(x) return(rbind(leg1(x),leg2(x),leg3(x),leg4(x),leg5(x),leg6(x),leg7(x),leg8(x),leg9(x),leg10(x))) #generate functional time series ( for simulations) rmuln<-function(d,x=1/2){ X=rnorm(d) D=diag(rep(1,d)) D[abs(row(D) - col(D)) == 1] <- x return(D%*%X) } rmult<-function(d,x=1/2,f=8){ X=rt(d,df=f)/sqrt(f/(f-2)) D=diag(rep(1,d)) D[abs(row(D) - col(D)) == 1] <- x return(D%*%X) } #generate functional innovations ( for simulations) gnfourier<-function(d,x=1/2,decay=(1:d)^(-2),res=100){ t=(1:res)/res X<-rmuln(d,x) return(as.vector(t(X*decay)%*%(cos((2*pi*(1:d))%*%t(t))+sin((2*pi*(1:d))%*%t(t))))) } gtfourier<-function(d,x=1/2,decay=(1:d)^(-2),res=100,f=8){ t=(1:res)/res X<-rmult(d,x,f) return(as.vector(t(X*decay)%*%(cos((2*pi*(1:d))%*%t(t))+sin((2*pi*(1:d))%*%t(t))))) } gnlegendre<-function(d,x=1/2,decay=(1:d)^(-2),res=100){ t=(1:res)/res X<-rmuln(d,x) return(as.vector(t(X*decay)%*%(leg(t)[1:d,]))) } gtlegendre<-function(d,x=1/2,decay=(1:d)^(-2),res=100,f=8){ t=(1:res)/res X<-rmult(d,x,f) return(as.vector(t(X*decay)%*%(leg(t)[1:d,]))) } #############simple functional linear time series models ( for simulations) gfma1n<-function(d,n,x=1/2,decay=(1:d)^(-2),res=100,a=1/2){ aa=matrix(0,nrow=(n+1),ncol=res) for(i in 1:(n+1)) aa[i,]=gnfourier(d,x,decay,res) return(aa[2:(n+1),]+a*aa[1:n,]) } gfar1n<-function(d,n,x=1/2,decay=(1:d)^(-2),res=100,a=1/2){ aa=matrix(0,nrow=(n+50),ncol=res) aa[1,]=gnfourier(d,x,decay,res) for(i in 2:(n+50)){ inno<-gnfourier(d,x,decay,res) aa[i,]=inno+a*aa[i-1,] } return(aa[51:(n+50),]) } gfma1t<-function(d,n,x=1/2,decay=(1:d)^(-2),res=100,a=1/2,f=8){ aa=matrix(0,nrow=(n+1),ncol=res) for(i in 1:(n+1)) aa[i,]=gtfourier(d,x,decay,res,f) return(aa[2:(n+1),]+a*aa[1:n,]) } gfar1t<-function(d,n,x=1/2,decay=(1:d)^(-2),res=100,a=1/2,f=8){ aa=matrix(0,nrow=(n+50),ncol=res) aa[1,]=gtfourier(d,x,decay,res,f) for(i in 2:(n+50)){ inno<-gtfourier(d,x,decay,res,f) aa[i,]=inno+a*aa[i-1,] } return(aa[51:(n+50),]) } #####generated by Legendre polynomials ( for simulations) gfma1n_leg<-function(d,n,x=1/2,decay=(1:d)^(-2),res=100,a=1/2){ aa=matrix(0,nrow=(n+1),ncol=res) for(i in 1:(n+1)) aa[i,]=gnlegendre(d,x,decay,res) return(aa[2:(n+1),]+a*aa[1:n,]) } gfar1n_leg<-function(d,n,x=1/2,decay=(1:d)^(-2),res=100,a=1/2){ aa=matrix(0,nrow=(n+50),ncol=res) aa[1,]=gnlegendre(d,x,decay,res) for(i in 2:(n+50)){ inno<-gnlegendre(d,x,decay,res) aa[i,]=inno+a*aa[i-1,] } return(aa[51:(n+50),]) } gfma1t_leg<-function(d,n,x=1/2,decay=(1:d)^(-2),res=100,a=1/2,f=8){ aa=matrix(0,nrow=(n+1),ncol=res) for(i in 1:(n+1)) aa[i,]=gtlegendre(d,x,decay,res,f) return(aa[2:(n+1),]+a*aa[1:n,]) } gfar1t_leg<-function(d,n,x=1/2,decay=(1:d)^(-2),res=100,a=1/2,f=8){ aa=matrix(0,nrow=(n+50),ncol=res) aa[1,]=gtlegendre(d,x,decay,res,f) for(i in 2:(n+50)){ inno<-gtlegendre(d,x,decay,res,f) aa[i,]=inno+a*aa[i-1,] } return(aa[51:(n+50),]) } ################generate simple panel functional time series ( for simulations) gpfar1n<-function(p,d,n,x=1/2,decay=(1:d)^(-2),res=100,a=1/2,pdep=1/2){ aa<-array(0,dim=c(p,n,res)) bb<-array(0,dim=c(p,n,res)) for(i in 1:p) aa[i,,]=gfar1n(d,n,x=1/2,decay,res,a) bb[1,,]=aa[1,,]+pdep*aa[2,,] for(i in 2:(p-1)) bb[i,,]=aa[i,,]+pdep*(aa[i-1,,]+aa[i+1,,]) bb[p,,]=aa[p,,]+pdep*aa[(p-1),,] rm(aa) return(bb) } gpfma1n<-function(p,d,n,x=1/2,decay=(1:d)^(-2),res=100,a=1/2,pdep=1/2){ aa<-array(0,dim=c(p,n,res)) bb<-array(0,dim=c(p,n,res)) for(i in 1:p) aa[i,,]=gfma1n(d,n,x=1/2,decay,res,a) bb[1,,]=aa[1,,]+pdep*aa[2,,] for(i in 2:(p-1)) bb[i,,]=aa[i,,]+pdep*(aa[i-1,,]+aa[i+1,,]) bb[p,,]=aa[p,,]+pdep*aa[(p-1),,] rm(aa) return(bb) } gpfar1t<-function(p,d,n,x=1/2,decay=(1:d)^(-2),res=100,a=1/2,pdep=1/2,f=8){ aa<-array(0,dim=c(p,n,res)) bb<-array(0,dim=c(p,n,res)) for(i in 1:p) aa[i,,]=gfar1t(d,n,x=1/2,decay,res,a,f) bb[1,,]=aa[1,,]+pdep*aa[2,,] for(i in 2:(p-1)) bb[i,,]=aa[i,,]+pdep*(aa[i-1,,]+aa[i+1,,]) bb[p,,]=aa[p,,]+pdep*aa[(p-1),,] rm(aa) return(bb) } gpfma1t<-function(p,d,n,x=1/2,decay=(1:d)^(-2),res=100,a=1/2,pdep=1/2,f=8){ aa<-array(0,dim=c(p,n,res)) bb<-array(0,dim=c(p,n,res)) for(i in 1:p) aa[i,,]=gfma1t(d,n,x=1/2,decay,res,a,f) bb[1,,]=aa[1,,]+pdep*aa[2,,] for(i in 2:(p-1)) bb[i,,]=aa[i,,]+pdep*(aa[i-1,,]+aa[i+1,,]) bb[p,,]=aa[p,,]+pdep*aa[(p-1),,] rm(aa) return(bb) } ###Legendre polynomials generated panel functional time series ( for simulations) gpfar1n_leg<-function(p,d,n,x=1/2,decay=(1:d)^(-2),res=100,a=1/2,pdep=1/2){ aa<-array(0,dim=c(p,n,res)) bb<-array(0,dim=c(p,n,res)) for(i in 1:p) aa[i,,]=gfar1n_leg(d,n,x=1/2,decay,res,a) bb[1,,]=aa[1,,]+pdep*aa[2,,] for(i in 2:(p-1)) bb[i,,]=aa[i,,]+pdep*(aa[i-1,,]+aa[i+1,,]) bb[p,,]=aa[p,,]+pdep*aa[(p-1),,] rm(aa) return(bb) } gpfma1n_leg<-function(p,d,n,x=1/2,decay=(1:d)^(-2),res=100,a=1/2,pdep=1/2){ aa<-array(0,dim=c(p,n,res)) bb<-array(0,dim=c(p,n,res)) for(i in 1:p) aa[i,,]=gfma1n_leg(d,n,x=1/2,decay,res,a) bb[1,,]=aa[1,,]+pdep*aa[2,,] for(i in 2:(p-1)) bb[i,,]=aa[i,,]+pdep*(aa[i-1,,]+aa[i+1,,]) bb[p,,]=aa[p,,]+pdep*aa[(p-1),,] rm(aa) return(bb) } gpfar1t_leg<-function(p,d,n,x=1/2,decay=(1:d)^(-2),res=100,a=1/2,pdep=1/2,f=8){ aa<-array(0,dim=c(p,n,res)) bb<-array(0,dim=c(p,n,res)) for(i in 1:p) aa[i,,]=gfar1t_leg(d,n,x=1/2,decay,res,a,f) bb[1,,]=aa[1,,]+pdep*aa[2,,] for(i in 2:(p-1)) bb[i,,]=aa[i,,]+pdep*(aa[i-1,,]+aa[i+1,,]) bb[p,,]=aa[p,,]+pdep*aa[(p-1),,] rm(aa) return(bb) } gpfma1t_leg<-function(p,d,n,x=1/2,decay=(1:d)^(-2),res=100,a=1/2,pdep=1/2,f=8){ aa<-array(0,dim=c(p,n,res)) bb<-array(0,dim=c(p,n,res)) for(i in 1:p) aa[i,,]=gfma1t_leg(d,n,x=1/2,decay,res,a,f) bb[1,,]=aa[1,,]+pdep*aa[2,,] for(i in 2:(p-1)) bb[i,,]=aa[i,,]+pdep*(aa[i-1,,]+aa[i+1,,]) bb[p,,]=aa[p,,]+pdep*aa[(p-1),,] rm(aa) return(bb) } ###Generate VAR(1)model ( for simulations) gpfvar1t_leg<-function(p,d,n,x=1/2,decay=(1:d)^(-2),res=100,a=1/2,f=8){ inno<-array(0,dim=c(p,n+50,res)) for(i in 1:p){ for(j in 1:(n+50)) inno[i,j,]=gtlegendre(d,x,decay,res,f) } D=diag(rep(1,p)) D[abs(row(D) - col(D)) == 1] <- 1/2 aa<-array(0,dim=c(p,n+50,res)) aa[,1,]<-inno[,1,] for(j in 2:(n+50)) aa[,j,]=a*(D%*%aa[,j-1,])+inno[,j,] return(aa[,51:(n+50),]) } ##########data analysis r=40 n=103 res=50 #dat is the data array dat<-array(0,dim=c(r,n,res)) dd<-read.csv(file = 'mean_Brampton_1879.csv') dat[1,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Brandon_1901.csv') dat[2,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Burnaby_1895.csv') dat[3,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Calgary_1885.csv') dat[4,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Charlottetown_1889.csv') dat[5,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Dawson_1901.csv') dat[6,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Edmonton_1882.csv') dat[7,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Fredericton_1907.csv') dat[8,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Halifax_1894.csv') dat[9,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Hamilton_1876.csv') dat[10,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Kamloops_1901.csv') dat[11,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Kelowna_1903.csv') dat[12,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Kenora_1916.csv') dat[13,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Kitchener_1874.csv') dat[14,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_La tuque_1904.csv') dat[15,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Laval_1874.csv') dat[16,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Lethbridge_1902.csv') dat[17,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_London_1868.csv') dat[18,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Markham_1897.csv') dat[19,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Medicine Hat_1912.csv') dat[20,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Mississauga_1841.csv') dat[21,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Montreal_1872.csv') dat[22,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Moose Jaw_1901.csv') dat[23,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Ottawa_1883.csv') dat[24,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Peterborough_1900.csv') dat[25,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Prince George_1916.csv') dat[26,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Prince Rupert_1911.csv') dat[27,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Quebec_1877.csv') dat[28,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Quesnel_1904.csv') dat[29,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Saint John_1900.csv') dat[30,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Saskatoon_1902.csv') dat[31,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Surrey_1895.csv') dat[32,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Sydney_1901.csv') dat[33,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Timmins_1914.csv') dat[34,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Toronto_1841.csv') dat[35,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Vancouver_1895.csv') dat[36,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Vaughan_1841.csv') dat[37,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Victoria_1891.csv') dat[38,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Windsor_1897.csv') dat[39,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) dd<-read.csv(file = 'mean_Yarmouth_1901.csv') dat[40,,]=t(dd[,(dim(dd)[2]-102):dim(dd)[2]]) #JSCB construction m=mv(dat,mrange=1:15) aaa=jscb(dat,m=10,repli=10000,alpha=c(0.001,0.01,0.05,0.1)) xx=1:50/50 yy1<-cos(2*pi*xx) zz1<-sin(2*pi*xx) yy2<-cos(4*pi*xx) zz2<-sin(4*pi*xx) #fitting sinusoidal trends par(mfrow=c(3,5)) fit1<-lm(aaa$fmean[4,]~yy1+zz1+yy2+zz2) plot(xx,fit1$fitted,lty=3,type="l",xlab="",ylab="mean_temp",ylim=c(-15,20),main="Calgary") lines(xx,aaa$upper[4,,3]) lines(xx,aaa$lower[4,,3]) fit1<-lm(aaa$fmean[6,]~yy1+zz1+yy2+zz2) plot(xx,fit1$fitted,lty=3,type="l",xlab="",ylab="mean_temp",ylim=c(-33,15),main="Dawson") lines(xx,aaa$upper[6,,3]) lines(xx,aaa$lower[6,,3]) fit1<-lm(aaa$fmean[7,]~yy1+zz1+yy2+zz2) plot(xx,fit1$fitted,lty=3,type="l",xlab="",ylab="mean_temp",ylim=c(-18,20),main="Edmonton") lines(xx,aaa$upper[7,,3]) lines(xx,aaa$lower[7,,3]) fit1<-lm(aaa$fmean[8,]~yy1+zz1+yy2+zz2) plot(xx,fit1$fitted,lty=3,type="l",xlab="",ylab="mean_temp",ylim=c(-13,20),main="Fredericton") lines(xx,aaa$upper[8,,3]) lines(xx,aaa$lower[8,,3]) fit1<-lm(aaa$fmean[9,]~yy1+zz1+yy2+zz2) plot(xx,fit1$fitted,lty=3,type="l",xlab="",ylab="mean_temp",ylim=c(-8,20),main="Halifax") lines(xx,aaa$upper[9,,3]) lines(xx,aaa$lower[9,,3]) fit1<-lm(aaa$fmean[16,]~yy1+zz1+yy2+zz2) plot(xx,fit1$fitted,lty=3,type="l",xlab="",ylab="mean_temp",ylim=c(-12,21),main="Laval") lines(xx,aaa$upper[16,,3]) lines(xx,aaa$lower[16,,3]) fit1<-lm(aaa$fmean[22,]~yy1+zz1+yy2+zz2) plot(xx,fit1$fitted,lty=3,type="l",xlab="",ylab="mean_temp",ylim=c(-14,23),main="Montreal") lines(xx,aaa$upper[22,,3]) lines(xx,aaa$lower[22,,3]) fit1<-lm(aaa$fmean[24,]~yy1+zz1+yy2+zz2) plot(xx,fit1$fitted,lty=3,type="l",xlab="",ylab="mean_temp",ylim=c(-14,23),main="Ottawa") lines(xx,aaa$upper[24,,3]) lines(xx,aaa$lower[24,,3]) fit1<-lm(aaa$fmean[27,]~yy1+zz1+yy2+zz2) plot(xx,fit1$fitted,lty=3,type="l",xlab="",ylab="mean_temp",ylim=c(-1,14),main="Prince Rupert") lines(xx,aaa$upper[27,,3]) lines(xx,aaa$lower[27,,3]) fit1<-lm(aaa$fmean[30,]~yy1+zz1+yy2+zz2) plot(xx,fit1$fitted,lty=3,type="l",xlab="",ylab="mean_temp",ylim=c(-10,17),main="Saint John") lines(xx,aaa$upper[30,,3]) lines(xx,aaa$lower[30,,3]) fit1<-lm(aaa$fmean[31,]~yy1+zz1+yy2+zz2) plot(xx,fit1$fitted,lty=3,type="l",xlab="",ylab="mean_temp",ylim=c(-22,20),main="Saskatoon") lines(xx,aaa$upper[31,,3]) lines(xx,aaa$lower[31,,3]) fit1<-lm(aaa$fmean[35,]~yy1+zz1+yy2+zz2) plot(xx,fit1$fitted,lty=3,type="l",xlab="",ylab="mean_temp",ylim=c(-8,23),main="Toronto") lines(xx,aaa$upper[35,,3]) lines(xx,aaa$lower[35,,3]) fit1<-lm(aaa$fmean[36,]~yy1+zz1+yy2+zz2) plot(xx,fit1$fitted,lty=3,type="l",xlab="",ylab="mean_temp",ylim=c(-1,20),main="Vancouver") lines(xx,aaa$upper[36,,3]) lines(xx,aaa$lower[36,,3]) fit1<-lm(aaa$fmean[38,]~yy1+zz1+yy2+zz2) plot(xx,fit1$fitted,lty=3,type="l",xlab="",ylab="mean_temp",ylim=c(2,17),main="Victoria") lines(xx,aaa$upper[38,,3]) lines(xx,aaa$lower[38,,3]) fit1<-lm(aaa$fmean[39,]~yy1+zz1+yy2+zz2) plot(xx,fit1$fitted,lty=3,type="l",xlab="",ylab="mean_temp",ylim=c(-8,23.5),main="Windsor") lines(xx,aaa$upper[39,,3]) lines(xx,aaa$lower[39,,3]) ############### test of parallelism aa<-Fptest(dat,m=10,repli=10000,alpha=c(0.001,0.002,0.005,0.008,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1)) ###calculating parallelism test statistic for each pair of locations r=dim(dat)[1] n=dim(dat)[2] fres=dim(dat)[3] cdat<-apply(dat,c(1,2),mean) transdat=dat for(i in 1:fres) transdat[,,i]=dat[,,i]-cdat rm(cdat) xx<-array(0,dim=c(r*(r-1),n,fres)) for(i in 1:r){ for(j in 1:r){ if(ji) xx[(r-1)*(i-1)+j-1,,]=transdat[i,,]-transdat[j,,] if(j==i) next } } rm(transdat) stdx1=std(xx) mm=apply(xx,c(1,3),mean) rm(xx) Tpair=sqrt(n)*apply(abs(mm)/stdx1,1,max) ###calculating parallelism test statistics for 15 selected cities a=10 Tpair[(a-1)*(r-1)+18-1] Tpair[(a-1)*(r-1)+14-1] Tpair[(a-1)*(r-1)+35-1] Tpair[(a-1)*(r-1)+39-1] Tpair[(a-1)*(r-1)+28-1] Tpair[(a-1)*(r-1)+24-1] Tpair[(a-1)*(r-1)+22-1] Tpair[(a-1)*(r-1)+4] Tpair[(a-1)*(r-1)+7] Tpair[(a-1)*(r-1)+5] Tpair[(a-1)*(r-1)+9] Tpair[(a-1)*(r-1)+38-1] Tpair[(a-1)*(r-1)+36-1] Tpair[(a-1)*(r-1)+31-1] a=18 Tpair[(a-1)*(r-1)+14] Tpair[(a-1)*(r-1)+35-1] Tpair[(a-1)*(r-1)+39-1] Tpair[(a-1)*(r-1)+28-1] Tpair[(a-1)*(r-1)+24-1] Tpair[(a-1)*(r-1)+22-1] Tpair[(a-1)*(r-1)+4] Tpair[(a-1)*(r-1)+7] Tpair[(a-1)*(r-1)+5] Tpair[(a-1)*(r-1)+9] Tpair[(a-1)*(r-1)+38-1] Tpair[(a-1)*(r-1)+36-1] Tpair[(a-1)*(r-1)+31-1] a=14 Tpair[(a-1)*(r-1)+35-1] Tpair[(a-1)*(r-1)+39-1] Tpair[(a-1)*(r-1)+28-1] Tpair[(a-1)*(r-1)+24-1] Tpair[(a-1)*(r-1)+22-1] Tpair[(a-1)*(r-1)+4] Tpair[(a-1)*(r-1)+7] Tpair[(a-1)*(r-1)+5] Tpair[(a-1)*(r-1)+9] Tpair[(a-1)*(r-1)+38-1] Tpair[(a-1)*(r-1)+36-1] Tpair[(a-1)*(r-1)+31-1] a=35 Tpair[(a-1)*(r-1)+39-1] Tpair[(a-1)*(r-1)+28] Tpair[(a-1)*(r-1)+24] Tpair[(a-1)*(r-1)+22] Tpair[(a-1)*(r-1)+4] Tpair[(a-1)*(r-1)+7] Tpair[(a-1)*(r-1)+5] Tpair[(a-1)*(r-1)+9] Tpair[(a-1)*(r-1)+38-1] Tpair[(a-1)*(r-1)+36-1] Tpair[(a-1)*(r-1)+31] a=39 Tpair[(a-1)*(r-1)+28] Tpair[(a-1)*(r-1)+24] Tpair[(a-1)*(r-1)+22] Tpair[(a-1)*(r-1)+4] Tpair[(a-1)*(r-1)+7] Tpair[(a-1)*(r-1)+5] Tpair[(a-1)*(r-1)+9] Tpair[(a-1)*(r-1)+38] Tpair[(a-1)*(r-1)+36] Tpair[(a-1)*(r-1)+31] a=28 Tpair[(a-1)*(r-1)+24] Tpair[(a-1)*(r-1)+22] Tpair[(a-1)*(r-1)+4] Tpair[(a-1)*(r-1)+7] Tpair[(a-1)*(r-1)+5] Tpair[(a-1)*(r-1)+9] Tpair[(a-1)*(r-1)+38-1] Tpair[(a-1)*(r-1)+36-1] Tpair[(a-1)*(r-1)+31-1] a=24 Tpair[(a-1)*(r-1)+22] Tpair[(a-1)*(r-1)+4] Tpair[(a-1)*(r-1)+7] Tpair[(a-1)*(r-1)+5] Tpair[(a-1)*(r-1)+9] Tpair[(a-1)*(r-1)+38-1] Tpair[(a-1)*(r-1)+36-1] Tpair[(a-1)*(r-1)+31-1] a=22 Tpair[(a-1)*(r-1)+4] Tpair[(a-1)*(r-1)+7] Tpair[(a-1)*(r-1)+5] Tpair[(a-1)*(r-1)+9] Tpair[(a-1)*(r-1)+38-1] Tpair[(a-1)*(r-1)+36-1] Tpair[(a-1)*(r-1)+31-1] a=4 Tpair[(a-1)*(r-1)+7-1] Tpair[(a-1)*(r-1)+5-1] Tpair[(a-1)*(r-1)+9-1] Tpair[(a-1)*(r-1)+38-1] Tpair[(a-1)*(r-1)+36-1] Tpair[(a-1)*(r-1)+31-1] a=7 Tpair[(a-1)*(r-1)+5] Tpair[(a-1)*(r-1)+9-1] Tpair[(a-1)*(r-1)+38-1] Tpair[(a-1)*(r-1)+36-1] Tpair[(a-1)*(r-1)+31-1] a=5 Tpair[(a-1)*(r-1)+9-1] Tpair[(a-1)*(r-1)+38-1] Tpair[(a-1)*(r-1)+36-1] Tpair[(a-1)*(r-1)+31-1] a=9 Tpair[(a-1)*(r-1)+38-1] Tpair[(a-1)*(r-1)+36-1] Tpair[(a-1)*(r-1)+31-1] a=38 Tpair[(a-1)*(r-1)+36] Tpair[(a-1)*(r-1)+31] a=36 Tpair[(a-1)*(r-1)+31] #######JSCB for Victoria-Vancouver and Calgary-Edmonton xx<-array(0,dim=c(2,n,fres)) xx[1,,]=dat[4,,]-dat[7,,] xx[2,,]=dat[38,,]-dat[36,,] hh<-jscb(xx,m=10,repli=10000,alpha=c(0.01,0.05)) xxx<-1:50/50 par(mfrow=c(1,2)) plot(xxx,hh$upper[1,,2],type="l",xlab="",ylab="temp.",ylim=c(-2.5,5),main="Victoria-Vancouver") lines(xxx,hh$lower[1,,2]) plot(xxx,hh$upper[2,,2],type="l",xlab="",ylab="temp.",ylim=c(-3,3),main="Calgary-Edmonton") lines(xxx,hh$lower[2,,2]) ###################ACF plots of 4 representative locations Tor=dat[35,,] Mon<-dat[22,,] Van<-dat[36,,] Cal<-dat[4,,] t1=15 t2=30 t3=45 par(mfrow=c(4,3)) xx=acf(Tor[,t1],plot=FALSE) plot(xx, main="Toronto, u=0.3") xx=acf(Tor[,t2],plot=FALSE) plot(xx, main="Toronto, u=0.6") xx=acf(Tor[,t3],plot=FALSE) plot(xx, main="Toronto, u=0.9") xx=acf(Mon[,t1],plot=FALSE) plot(xx, main="Montreal, u=0.3") xx=acf(Mon[,t2],plot=FALSE) plot(xx, main="Montreal, u=0.6") xx=acf(Mon[,t3],plot=FALSE) plot(xx, main="Montreal, u=0.9") xx=acf(Van[,t1],plot=FALSE) plot(xx, main="Vancouver, u=0.3") xx=acf(Van[,t2],plot=FALSE) plot(xx, main="Vancouver, u=0.6") xx=acf(Van[,t3],plot=FALSE) plot(xx, main="Vancouver, u=0.9") xx=acf(Cal[,t1],plot=FALSE) plot(xx, main="Calgary, u=0.3") xx=acf(Cal[,t2],plot=FALSE) plot(xx, main="Calgary, u=0.6") xx=acf(Cal[,t3],plot=FALSE) plot(xx, main="Calgary, u=0.9")