## Each variable should be TxN ## The regressor matrix is saved as a 3D array TxNxK ## Strictly exogenous variables are processed by LSDV ## Weakly exogenous variables are processed by FOD/BOD or FOD/x[it] ## rho is estimated by HP ## ## Z <- array(1:30,c(3,5,2)) # TxNxK ## z <- matrix(z,3*5,2) # NTxK ## ## vec is an NT-vector ## mat is TxN ## dat is NTxK ## arr is TxNxK DIM <- function(x) { if (is.vector(x)) length(x) else dim(x) } ## Functions to process Stata-type data get.t.set <- function(ivec,tvec) { x <- NULL for (i in unique(ivec)) { ti <- tvec[i==ivec] x <- union(x,ti) } x } get.std.tvec <- function(tvec, t.set) { x <- rep(NA,length(tvec)) s <- 1 for (t in t.set) { x[t==tvec] <- s s <- s+1 } x } convert.vec.to.mat <- function(x,ivec,tvec,t.set=NULL) { # to NxT if (is.null(t.set)) t.set <- get.t.set(ivec,tvec) nn <- length(unique(ivec)) tt <- length(t.set) stdt <- get.std.tvec(tvec,t.set) z <- matrix(NA,tt,nn) for (i in unique(ivec)) { ii <- i==ivec ti <- tvec[ii] zi <- rep(NA,tt) zi[stdt] <- x[ii] z[,i] <- zi } z } to.regular.panel <- function(x,ivec,tvec,t.set=NULL) { if (is.vector(x)) { to.vec(convert.vec.to.mat(x,ivec,tvec,t.set)) } else if (is.matrix(x)) { z <- matrix(NA,NROW(x),NCOL(x)) for (j in 1:NCOL(x)) z[,j] <- to.vec(convert.vec.to.mat(x[,j],ivec,tvec,t.set)) z } } convert.dat.to.arr <- function(x,ivec,tvec,t.set=NULL) { # dat to arr if (is.vector(x)) { read.panel.vec(x,ivec,tvec,t.set) } else { if (is.null(t.set)) t.set <- get.t.set(ivec,tvec) nn <- length(unique(ivec)) tt <- length(t.set) stdt <- get.std.tvec(tvec,t.set) } } is.arr <- function(x) length(dim(x))==3 ## to.vec works for mat, dat and arr to.vec <- function(x) as.vector(x) vec.to.mat <- function(x,n,t=round(length(x)/n)) matrix(x,t,n) ## vec.to.arr works for vec vec.to.arr <- function(x,n,k,t=round(length(x)/n/k)) { array(x,c(t,n,k)) } ## mat.to.arr converts mat to arr [,,1] mat.to.arr <- function(x) { dim(x) <- c(NROW(x),NCOL(x),1) x } dat.to.arr <- function(x,n,t=round(NROW(x)/n)) { z <- array(NA,c(t,n,NCOL(x))) x <- matrix(x,NROW(x),NCOL(x)) for (k in 1:NCOL(x)) z[,,k] <- vec.to.mat(x[,k],n,t) z } arr.to.dat <- function(x) { d <- dim(x) if (length(d)<3) to.vec(x) else matrix(x,d[1]*d[2],d[3]) } mat.within <- function(x) { x-rep(1,NROW(x))%o%colMeans(x) } arr.within <- function(x) { if (!is.3d(x)) x <- array(x,c(NROW(x),NCOL(x),1)) d <- dim(x) z <- array(0,d) for (k in 1:(d[3])) z[,,k] <- mat.within(x[,,k]) z } vec.within <- function(x,n) to.vec(mat.within(vec.to.mat(x,n))) dat.within <- function(x,n) arr.to.dat(arr.within(dat.to.arr(x,n))) mat.trimr <- function(x,top,bot) x[(1+top):(NROW(x)-bot),] arr.trimr <- function(x,top,bot) { if (top==0 && bot==0) { x } else { x[(1+top):(dim(x)[1]-bot),,] } } vec.trimr <- function(x,n,top,bot) { to.vec(mat.trimr(vec.to.mat(x,n),top,bot)) } dat.trimr <- function(x,n,top,bot) { arr.to.dat(arr.trimr(dat.to.arr(x,n),top,bot)) } mat.diff <- function(x,rho=1) mat.trimr(x,1,0)-rho*mat.trimr(x,0,1) arr.diff <- function(x,rho=1) arr.trimr(x,1,0)-rho*arr.trimr(x,0,1) vec.diff <- function(x,n,k=NULL,rho=1) { if (is.null(k)) { # vec to.vec(mat.diff(vec.to.mat(x,n),rho=rho)) } else { to.vec(arr.diff(vec.to.arr(x,n,k),rho=rho)) } } dat.diff <- function(x,n,rho=1) { arr.to.dat(arr.diff(dat.to.arr(x,n),rho=rho)) } arr.lag <- function(x,lag=1) arr.trimr(x,0,1) panel.default.lags <- function(x) { n <- NCOL(x)-sum(colSums(is.na(x))==NROW(x)) 1:n-1 } get.fod <- function(n) { # XF' is FOD transformation if (n<2) { NULL } else { x <- matrix(0,n-1,n) for (i in 1:(n-1)) { x[i,i] <- 1 x[i,(i+1):n] <- -1/(n-i) } sqrt((n:2-1)/n:2)*x } } get.bod <- function(n) { # XB' is BOD transformation if (n<2) { NULL } else { x <- matrix(0,n-1,n) for (i in 1:(n-1)) { x[i,i+1] <- 1 x[i,1:i] <- -1/i } x } } arr.filter <- function(x,y) { d <- dim(x) d[1] <- NROW(y) z <- array(0,d) for (k in 1:(d[3])) z[,,k] <- y%*%x[,,k] z } arr.fod <- function(x,fod=NULL) { if (is.null(fod)) fod <- get.fod(NROW(x[,,1])) arr.filter(x,fod) } arr.bod <- function(x,bod=NULL) { if (is.null(bod)) bod <- get.bod(NROW(x[,,1])) arr.filter(x,bod) } dat.fod <- function(x,n,fod=NULL) { arr.to.dat(arr.fod(dat.to.arr(x,n),fod)) } dat.bod <- function(x,n,bod=NULL) { arr.to.dat(arr.bod(dat.to.arr(x,n),bod)) } is.3d <- function(x) length(dim(x))==3 panel.fod <- function(x,n=NULL,fod=NULL,k=NULL) { if (is.vector(x)) { to.vec(arr.fod(vec.to.arr(x,n,k), fod)) } else if (is.3d(x)) { arr.fod(x,fod) } else if (is.null(n)) { # mat if (is.null(fod)) fod <- get.fod(NROW(x)) fod%*%x } else { # dat dat.fod(x,n,fod) } } panel.bod <- function(x,bod=NULL,n=NULL,k=NULL) { if (is.vector(x)) { to.vec(arr.bod(vec.to.arr(x,n,k), bod)) } else if (is.3d(x)) { arr.bod(x,bod) } else if (is.null(n)) { # mat if (is.null(bod)) bod <- get.bod(NROW(x)) bod%*%x } else { # dat dat.bod(x,n,bod) } } .xthp.nox <- function(dy,n=NULL) { # dy is NTx1 or TxN if (!is.null(n)) dy <- vec.to.mat(dy,n) x <- mat.trimr(dy,0,1) z <- 2*mat.trimr(dy,1,0)+x rho <- sum(x*z)/sum(x^2) names(rho) <- NULL rho } .xthp.nox.all <- function(dy,n=NULL,rho=NULL,x=NULL,z=NULL,sxx=NULL) { if (!is.null(n)) dy <- vec.to.mat(dy,n) if (is.null(x)) x <- mat.trimr(dy,0,1) if (is.null(z)) z <- 2*mat.trimr(dy,1,0)+x if (is.null(sxx)) sxx <- sum(x^2) if (is.null(rho)) rho <- sum(x*z)/sxx n <- NCOL(z) t <- NROW(z) ee <- z-rho*x v2 <- 2*(1+rho)/(n*t) v1 <- crossprod(colSums(x*ee))/(sxx^2) #*((t+2)/t) se <- sqrt(v1) t.rho <- rho/se p.rho <- 2*pnorm(-abs(t.rho)) ans <- matrix(c(rho, se, t.rho, p.rho),1,4) rownames(ans) <- "L1" colnames(ans) <- c("Coef.", "Std. Err.", "z", "P>|z|") list(rho=rho, se=sqrt(c(v1,v2)), n=n, df=n*t, coefficients=ans) } .force.fmt <- function(x,n=NULL) { if (is.vector(x)) { x <- if (!is.null(n)) matrix(x,NROW(x),NCOL(x)) else vec.to.arr(x,n,1) } x } .xthp.x <- function(dy,dx,yd2,xd2,yd1,xd1,n=NULL,iter=25,eps=1e-6) { ## if n==NULL, then dy:mat, dx:arr, yd:mat, xd:arr ## if n!=NULL, then dy:vec, dx:dat, yd:vec, xd:dat ## initial dx <- .force.fmt(dx,n) xd2 <- .force.fmt(xd2,n) xd1 <- .force.fmt(xd1,n) if (is.null(n)) { dy <- to.vec(dy) dx <- arr.to.dat(dx) yd1 <- to.vec(yd1) xd1 <- arr.to.dat(xd1) yd2 <- to.vec(yd2) xd2 <- arr.to.dat(xd2) } if (!is.matrix(xd1)) { xd1 <- matrix(xd1, NROW(xd1), NCOL(xd1)) xd2 <- matrix(xd2, NROW(xd2), NCOL(xd2)) } if (!is.matrix(dx)) dx <- matrix(dx, NROW(dx), NCOL(dx)) rho1 <- 0 rho0 <- rho1-Inf i <- 2 while (abs(rho1-rho0)>eps && i<=iter) { rho0 <- rho1 yy <- yd2-rho0*yd1 X <- xd2-rho0*xd1 lsdv <- lm(yy~X-1) b <- lsdv$coef e <- dy-as.vector(dx%*%b) rho1 <- .xthp.nox(e,n) i <- i+1 } nox <- .xthp.nox.all(e,n) sum.lsdv <- summary(lsdv) sigma <- sum.lsdv$sigma*sqrt(sum(sum.lsdv$df[1:2])/nox$df) cov.matrix <- sigma^2*sum.lsdv$cov.unscaled se.b <- sqrt(diag(cov.matrix)) t.val <- b/se.b p.val <- round(2*pnorm(-abs(t.val)),7) ans <- cbind(b,se.b,t.val,p.val) t.rho <- nox$rho/nox$se[1] p.rho <- round(2*pnorm(-abs(t.rho)),7) ans <- rbind(c(nox$rho, nox$se[1], t.rho, p.rho), ans) rownames(ans)[1] <- "L1" colnames(ans) <- c("Coef.", "Std. Err.", "z", "P>|z|") if (iter==1) rho0 <- NA list(rho=nox$rho, se.rho=nox$se, bhat=b, cov.bhat=cov.matrix, sigma.bhat=sigma, coefficients=ans, iter=i-1, eps=abs(rho1-rho0)) } .xthp.w <- function(dy,dw,yd2,wd2,yd1,wd1,zw,n=NULL,iter=25,eps=1e-6) { ## if n==NULL, then dy:mat, dx:arr, yd:mat, xd:arr ## if n!=NULL, then dy:vec, dx:dat, yd:vec, xd:dat ## initial dw <- .force.fmt(dw,n) wd2 <- .force.fmt(wd2,n) wd1 <- .force.fmt(wd1,n) if (is.null(n)) { dy <- to.vec(dy) dw <- arr.to.dat(dw) yd1 <- to.vec(yd1) wd1 <- arr.to.dat(wd1) yd2 <- to.vec(yd2) wd2 <- arr.to.dat(wd2) } if (!is.matrix(wd1)) { wd1 <- matrix(wd1, NROW(wd1), NCOL(wd1)) wd2 <- matrix(wd2, NROW(wd2), NCOL(wd2)) } if (!is.matrix(dw)) dw <- matrix(dw, NROW(dw), NCOL(dw)) rho1 <- 0 rho0 <- rho1-Inf i <- 2 while (abs(rho1-rho0)>eps && i<=iter) { rho0 <- rho1 yy <- yd2-rho0*yd1 W <- wd2-rho0*wd1 W <- lm(W~zw-1)$fitted ##yy <- lm(yy~zw-1)$fitted reg <- lm(yy~W-1) b <- reg$coef e <- dy-as.vector(dw%*%b) rho1 <- .xthp.nox(e,n) i <- i+1 } nox <- .xthp.nox.all(e,n) sum.reg <- summary(reg) yy <- yd2-rho1*yd1 W <- wd2-rho1*wd1 res <- yy-as.vector(W%*%b) sigma <- sqrt(sum(res^2)/sum(sum.reg$df[1:2])) cov.matrix <- sigma^2*sum.reg$cov.unscaled se.b <- sqrt(diag(cov.matrix)) t.val <- b/se.b p.val <- round(2*pnorm(-abs(t.val)),7) ans <- cbind(b,se.b,t.val,p.val) t.rho <- nox$rho/nox$se[1] p.rho <- round(2*pnorm(-abs(t.rho)),7) ans <- rbind(c(nox$rho, nox$se[1], t.rho, p.rho), ans) rownames(ans)[1] <- "L1" colnames(ans) <- c("Coef.", "Std. Err.", "z", "P>|z|") if (iter==1) rho0 <- NA list(rho=nox$rho, se.rho=nox$se, bhat=b, cov.bhat=cov.matrix, sigma.bhat=sigma, coefficients=ans, iter=i-1, eps=abs(rho1-rho0)) } .xthp.xw <- function(dy,dx,dw,yd2,xd2,wd2,yd1,xd1,wd1,zw, n=NULL,iter=25,eps=1e-5) { ## if n==NULL, then dy:mat, dx:arr, yd:mat, xd:arr ## if n!=NULL, then dy:vec, dx:dat, yd:vec, xd:dat ## initial dx <- .force.fmt(dx,n) xd2 <- .force.fmt(xd2,n) xd1 <- .force.fmt(xd1,n) wd2 <- .force.fmt(wd2,n) wd1 <- .force.fmt(wd1,n) if (is.null(n)) { dy <- to.vec(dy) dx <- arr.to.dat(dx) xd1 <- arr.to.dat(xd1) xd2 <- arr.to.dat(xd2) dw <- arr.to.dat(dw) yd1 <- to.vec(yd1) wd1 <- arr.to.dat(wd1) yd2 <- to.vec(yd2) wd2 <- arr.to.dat(wd2) } if (!is.matrix(xd1)) { xd1 <- matrix(xd1, NROW(xd1), NCOL(xd1)) xd2 <- matrix(xd2, NROW(xd2), NCOL(xd2)) } if (!is.matrix(wd1)) { wd1 <- matrix(wd1, NROW(wd1), NCOL(wd1)) wd2 <- matrix(wd2, NROW(wd2), NCOL(wd2)) } if (!is.matrix(dx)) dx <- matrix(dx, NROW(dx), NCOL(dx)) if (!is.matrix(dw)) dw <- matrix(dw, NROW(dw), NCOL(dw)) rho1 <- 0 rho0 <- rho1-Inf i <- 2 while (abs(rho1-rho0)>eps && i<=iter) { rho0 <- rho1 yy <- yd2-rho0*yd1 X <- xd2-rho0*xd1 W <- wd2-rho0*wd1 W <- lm(W~X+zw-1)$fitted ##yy <- lm(yy~X+zw-1)$fitted reg <- lm(yy~X+W-1) b <- reg$coef e <- dy-as.vector(cbind(dx,dw)%*%b) rho1 <- .xthp.nox(e,n) i <- i+1 } nox <- .xthp.nox.all(e,n) sum.reg <- summary(reg) yy <- yd2-rho1*yd1 X <- xd2-rho1*xd1 W <- wd2-rho1*wd1 res <- yy-as.vector(cbind(X,W)%*%b) sigma <- sqrt(sum(res^2)/sum(sum.reg$df[1:2])) cov.matrix <- sigma^2*sum.reg$cov.unscaled se.b <- sqrt(diag(cov.matrix)) t.val <- b/se.b p.val <- round(2*pnorm(-abs(t.val)),7) ans <- cbind(b,se.b,t.val,p.val) t.rho <- nox$rho/nox$se[1] p.rho <- round(2*pnorm(-abs(t.rho)),7) ans <- rbind(c(nox$rho, nox$se[1], t.rho, p.rho), ans) rownames(ans)[1] <- "L1" colnames(ans) <- c("Coef.", "Std. Err.", "z", "P>|z|") if (iter==1) rho0 <- NA list(rho=nox$rho, se.rho=nox$se, bhat=b, cov.bhat=cov.matrix, sigma.bhat=sigma, coefficients=ans, iter=i-1, eps=abs(rho1-rho0)) } xthp <- function(y,x=NULL,w=NULL,n,bod.iv=TRUE,iter=25,eps=1e-6) { ## data: Stata style ## y: depvar, x: exogenous, w: weakly exogenous y <- vec.to.mat(y,n) dy <- mat.diff(y) if (!is.null(x)) { x <- dat.to.arr(x,n) dx <- arr.to.dat(arr.diff(x)) } else { dx <- NULL } if (!is.null(w)) { w <- dat.to.arr(w,n) dw <- arr.to.dat(arr.diff(w)) } if (is.null(x) && is.null(w)) { .xthp.nox.all(dy) } else if (is.null(w)) { yd1 <- to.vec(mat.within(mat.trimr(y,0,1))) yd2 <- to.vec(mat.within(mat.trimr(y,1,0))) xd1 <- arr.to.dat(arr.within(arr.trimr(x,0,1))) xd2 <- arr.to.dat(arr.within(arr.trimr(x,1,0))) .xthp.x(dy,dx,yd2,xd2,yd1,xd1,n=n,iter=iter,eps=eps) } else { if (bod.iv) { fod <- get.fod(NROW(y)-1) bod <- get.bod(NROW(y)-1) yd1 <- to.vec(panel.fod(mat.trimr(y,0,1),fod=fod)) yd2 <- to.vec(panel.fod(mat.trimr(y,1,0),fod=fod)) wd1 <- arr.to.dat(panel.fod(arr.trimr(w,0,1),fod=fod)) wd2 <- arr.to.dat(panel.fod(arr.trimr(w,1,0),fod=fod)) zw <- arr.to.dat(panel.bod(arr.trimr(w,0,1),bod=bod)) } else { fod <- get.fod(NROW(y)-1) bod <- get.bod(NROW(y)-1) fod <- get.fod(NROW(y)-1) bod <- get.bod(NROW(y)-1) yd1 <- to.vec(panel.fod(mat.trimr(y,0,1),fod=fod)) yd2 <- to.vec(panel.fod(mat.trimr(y,1,0),fod=fod)) wd1 <- arr.to.dat(panel.fod(arr.trimr(w,0,1),fod=fod)) wd2 <- arr.to.dat(panel.fod(arr.trimr(w,1,0),fod=fod)) zw <- arr.to.dat(arr.trimr(w,0,2)) } if (is.null(x)) { .xthp.w(dy,dw,yd2,wd2,yd1,wd1,zw,n=n,iter=iter,eps=eps) } else { if (bod.iv) { xd1 <- arr.to.dat(panel.fod(arr.trimr(x,0,1),fod=fod)) xd2 <- arr.to.dat(panel.fod(arr.trimr(x,1,0),fod=fod)) } else { xd1 <- arr.to.dat(panel.fod(arr.trimr(x,0,1),fod=fod)) xd2 <- arr.to.dat(panel.fod(arr.trimr(x,1,0),fod=fod)) } .xthp.xw(dy,dx,dw,yd2,xd2,wd2,yd1,xd1,wd1,zw,n=n,iter=iter,eps=eps) } } }