| # ============================================================= | | # [FUNCTION]: OLS estimation for GSTAR(p;L1,…,Lp) models | | # ============================================================= | | # 3 dimension zeros matrix | | # –––––––––––––––- | | zeros <-function(m,n,p){ | | W<-rep(0, m*n*p) | | dim(W)<-c(m, n, p) | | W} | | # “vec” operator | | vec<-function(X){ | | a<-dim(X) | | Y<- t(X[1, ]) | | for (i in 2:a[1]){ | | Y<- cbind(Y,t(X[i, ]))} | | t(Y)} | | #–––––––––– | | # Inverse of matrix | | #–––––––––– | | inv<-function(X){ | | if(dim(X)[1]! = dim(X)[2]) stop(“THE MATRIX MUST BE SYMMETRIC!!!”) | | else{ | | if (det(X)==0) stop("THE MATRIX IS SINGULAR!!!") | | else { n<-dim(X)[1] | | solve(X) | | # ––––––––––––––––––––––––––––––– | | # Construction of vector Zi, for each i=1,…,N | | # ––––––––––––––––––––––––––– | | # Construction of vector Zi, for each i=1,…,N | | # ––––––––––––––––––––––––––– | | # suppose x = c(p,L1,..., Lp) represent the model order | | Zi<-function(Zt, x){ | | N<-dim(Zt)[1] #number of sites | | T<-dim(Zt)[2]-1 #number of time periods | | p<-x [1] | | Zi<-matrix(0, T-p+1, N) | | for (i in 1: N) | | Zi[,i]<-Zt[i,(p+1):(T+1)] | | Zi} | | # –––––––––––––––––––––––––––––––- | | # Construction of matrix Xi, for each i = 1,…,N | | # ––––––––––––––––––––––- | | Xi<-function(Zt,x) | | {N<-dim(Zt)[1] #number of sites | | T<-dim(Zt)[2]-1 #number of time periods | | p<-x[1] | | La<-x[2: length(x)] | | r<-lmd+1 # where lmd = the greatest order for weight matrices | | WZ<-zeros(N, T, r) | | for (k in 1: r) | | WZ[,, k]<-W [,, k]%*%Zt[, 1: T] | | Xi<-zeros((T-p+1),sum(La+1), N) | | if (p==1) | | { for (i in 1:N){ | | TR<-WZ[i,p:T,1:(La[1]+1)] | | Xi[,, i]<-TR | | if (p>=2){ | | for (i in 1:N){ | | TR<-WZ[i,p:T,1:(La[1]+1)] | | for (s in 2:p) | | TR<-cbind(TR,WZ[i,(p-s+1):(T-s+1),1:(La[s]+1)]) | | Xi[,, i]<-TR | | Xi} | | # –––––––––––––––––––––––––––––––- | | # OLS parameter of GSTAR model | | # –––––––––––––– | | gstar<-function(Zt,x){ | | p<-x[1] | | La<-x[2:length(x)] | | r<-lmd+1 # where lmd = the greatest order for weight matrices | | N<-dim(Zt)[1] #number of sites | | T<-dim(Zt)[2]-1 #number of time periods | | Xi<-Xi(Zt,x) | | Zi<-Zi(Zt,x) | | coef.OLS<-matrix(0,sum(La+1),N) | | col.name<-array(0,N) | | for (i in 1:N){ | | coef.OLS[, i]<-inv(t(Xi[,, i])%*%Xi[,,i])%*%t(Xi[,, i])%*%Zi[,i] | | col.name[i]<-paste("site",i)} | | colnames(coef.OLS)<-col.name | | round(coef.OLS,4)} | | # –––––––––––––––––––––––––––––––- | | # Residuals of GSTAR model | | # –––––––––––––– | | # (1). To find the LS estimates only, for example GSTAR(2;1,1), use | | # the command: | | # > gstar(Zt,c(2,1,1)) | | # where Zt is data matrix. | | # (2). To find the estimates, prediction values, and residuals | | # vector respectively, call the function by the following | | # commands: | | # > as.2<-res(Zt,c(2,1,1)) | | # > as.2$coef | | # > as.2$pred | | # > as.2$res | | # –––––––––––––––––––––––––––––––- | | res<-function(Zt,x){ | | coef<-gstar(Zt,x) | | Xi<-Xi(Zt,x) | | p<-x[1] | | La<-x[2:length(x)] | | N<-dim(Zt)[1] #number of sites | | T<-dim(Zt)[2]-1 #number of time periods | | Z.OLS<-matrix(0,T-p+1,N) | | res.OLS<-matrix(0,N,T-p+1) | | if (p==1){ | | for (i in 1:N){ | | if (La[1]!=0)Z.OLS[,i]<-Xi[,,i]%*%coef[,i] | | else Z.OLS[,i]<-Xi[,,i]*coef[,i] | | res.OLS[i,]<-t(Z.OLS[,i]-Zt[i,(p+1):(T+1)])} | |
}
| | if (p!=1){ | | for (i in 1:N){ | | Z.OLS[,i]<-Xi[,,i]%*%coef[,i] | | res.OLS[i,]<-t(Z.OLS[,i]-Zt[i,(p+1):(T+1)]) | | az<-new.env() | | az$Xi<-Xi # matrix Xi | | az$coef<-coef | | az$pred<-t(Z.OLS) | | az$res<-res.OLS | | ax<-as.list(az)} | | # –––––––––––––––––––––––––––––––- |
|
