[Splm-commits] r180 - pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Nov 23 23:56:34 CET 2013
Author: gpiras
Date: 2013-11-23 23:56:34 +0100 (Sat, 23 Nov 2013)
New Revision: 180
Modified:
pkg/R/ivplm.w2sls.R
pkg/R/likelihoodsFE.R
pkg/R/spfeml.R
pkg/R/spgm.R
pkg/R/spml.R
Log:
Lee and Yu with twoways
Modified: pkg/R/ivplm.w2sls.R
===================================================================
--- pkg/R/ivplm.w2sls.R 2013-11-06 23:12:32 UTC (rev 179)
+++ pkg/R/ivplm.w2sls.R 2013-11-23 22:56:34 UTC (rev 180)
@@ -121,7 +121,7 @@
endogwithin <-cbind(endogwithin, wywithin)
colnames(endogwithin)<-c(colnames(endog), "lambda")
-colnames(Xwithin)<-colnames(X)[-del]
+# colnames(Xwithin)<-colnames(X)[-del]
res<-spgm.tsls(ywithin, endogwithin, Xwithin, Hwithin)
Modified: pkg/R/likelihoodsFE.R
===================================================================
--- pkg/R/likelihoodsFE.R 2013-11-06 23:12:32 UTC (rev 179)
+++ pkg/R/likelihoodsFE.R 2013-11-23 22:56:34 UTC (rev 180)
@@ -22,7 +22,7 @@
}
-splaglm<-function(env, zero.policy = zero.policy, interval = interval, con = con, llprof = llprof, tol.solve= tol.solve, Hess = Hess, method = method){
+splaglm<-function(env, zero.policy = zero.policy, interval = interval, con = con, llprof = llprof, tol.solve= tol.solve, Hess = Hess, method = method, LeeYu = LeeYu, effects = effects){
xt <- get("xt", envir = env)
yt <- get("yt", envir = env)
@@ -30,6 +30,7 @@
con<-get("con", envir = env)
NT<-get("NT", envir = env)
T<-get("T", envir = env)
+n <- NT/T
listw<-get("listw", envir = env)
inde<-get("inde", envir = env)
interval1 <- get("interval1", envir = env)
@@ -49,8 +50,6 @@
opt <- optimize(conclikpan, interval = interval1, maximum = TRUE, env = env, tol = con$tol.opt)
-
-
#opt <- nlminb(0.02138744, conclikpan, lower = interval[1], upper= interval[2], env = env)
lambda <- opt$maximum
@@ -62,34 +61,74 @@
LL <- opt$objective
optres <- opt
- # nm <- paste(method, "opt", sep = "_")
- # timings[[nm]] <- proc.time() - .ptime_start
- # .ptime_start <- proc.time()
-
lm.lag <- lm((yt - lambda * wyt) ~ xt - 1)
p <- lm.lag$rank
r <- residuals(lm.lag)
fit <- yt - r
names(r) <- names(fit)
+ SSE <- crossprod(residuals(lm.lag))
+ s2 <- as.numeric(SSE)/NT
+
+if(LeeYu && effects == "spfe") s2 <- (T/(T-1)) * as.numeric(s2)
+if(LeeYu && effects == "tpfe") s2 <- (n/(n-1)) * as.numeric(s2)
+
betas <- coefficients(lm.lag)
names(betas) <- colnames(xt)
coefs <- c(lambda, betas)
- SSE <- deviance(lm.lag)
- s2 <- SSE/NT
- # coefsl <- c(s2, rho, betas)
-
- # timings[["coefs"]] <- proc.time() - .ptime_start
- # .ptime_start <- proc.time()
- # assign("first_time", TRUE, envir = env)
+
+if(LeeYu && effects == "sptpfe"){
+
+ tr <- function(A) sum(diag(A))
+ W <-listw2dgCMatrix(listw, zero.policy = zero.policy)
+ A <- solve(sparseMatrix(i=1:(NT/T), j=1:(NT/T), x=1) - lambda * W)
+ WA <- W %*% A
+ lag <- function(q) trash<-unlist(tapply(q,inde,function(TT) as.matrix(WA %*% TT), simplify=TRUE))
+ lag2 <- function(q) trash<-unlist(tapply(q,inde,function(TT) as.matrix(t(WA)%*%TT), simplify=TRUE))
+ WAxt <- apply(as.matrix(xt),2,lag)
+ WAWAxt<-apply(WAxt,2,lag2)
+ xtWAWAxt <- crossprod(xt,WAWAxt)
+ xtWAxt <- crossprod(xt,WAxt)
+ xtxt <- crossprod(xt)
+ one <- T*(tr(WA %*% WA) + tr(t(WA) %*% WA))
+ two <- 1/as.numeric(s2) * t(betas) %*% xtWAWAxt %*% betas
+ V <- one + two
+ zero <- rbind(rep(0, length(betas)))
+ col1 <- rbind(NT/(2 * (s2^2)), T*tr(WA)/s2, t(zero))
+ three <- (1/as.numeric(s2)) * xtWAxt %*% betas
+ col2 <- rbind(T*tr(WA)/s2, V, three )
+ col3 <- rbind(zero, t(three), 1/as.numeric(s2)* xtxt)
+ asyvar <- cbind(col1, col2, col3)
+ asyv <- solve(asyvar, tol = con$tol.solve)
+ rownames(asyv) <- colnames(asyv) <- c("sigma","lambda", colnames(xt))
+
+ init <- c((T/(T+1)), rep(1,p+1))
+
+ a3 <- rep(0,p)
+ a2 <- 1/(1 - lambda)
+ a1 <- 1/(2*s2)
+ a <- c(a1,a2,a3)
+ Bhat <- - (asyv/n) %*% a
+
+
+ coefs1 <- c(s2, lambda, betas)
+ Theta2 <- init * coefs1 + Bhat
+ betas <- as.numeric(Theta2[(3:(2+p))])
+ names(betas) <- colnames(xt)
+ lambda <-Theta2[2]
+ names(lambda) <- "lambda"
+ s2 <- Theta2[1]
+ coefs <- c(lambda, betas)
+
+}
### add numerical hessian
if(Hess){
- fd <- fdHess(coefs, f_sarpanel_hess, env)
+ fd <- fdHess(coefs, f_sarpanel_hess, env, LeeYu = LeeYu, effects = effects)
mat <- fd$Hessian
- fdHess<- solve(-(mat), tol.solve = tol.solve)
+ fdHess<- solve(-(mat), tol.solve = tol.solve)
rownames(fdHess) <- colnames(fdHess) <- c("lambda", colnames(xt))
lambda.se <- fdHess[1, 1]
@@ -105,20 +144,33 @@
W <-listw2dgCMatrix(listw, zero.policy = zero.policy)
A <- solve(sparseMatrix(i=1:(NT/T), j=1:(NT/T), x=1) - lambda * W)
WA <- W %*% A
- one <- T*(tr(WA %*% WA) + tr(t(WA) %*% WA))
-
- lag<-function(q) trash<-unlist(tapply(q,inde,function(TT) as.matrix(WA %*% TT), simplify=TRUE))
- lag2<-function(q) trash<-unlist(tapply(q,inde,function(TT) as.matrix(t(WA)%*%TT), simplify=TRUE))
- WAxt<-apply(as.matrix(xt),2,lag)
-# print(WAxt)
+ lag <- function(q) trash<-unlist(tapply(q,inde,function(TT) as.matrix(WA %*% TT), simplify=TRUE))
+ lag2 <- function(q) trash<-unlist(tapply(q,inde,function(TT) as.matrix(t(WA)%*%TT), simplify=TRUE))
+ WAxt <- apply(as.matrix(xt),2,lag)
WAWAxt<-apply(WAxt,2,lag2)
- # print(WAWAxt)
xtWAWAxt <- crossprod(xt,WAWAxt)
xtWAxt <- crossprod(xt,WAxt)
xtxt <- crossprod(xt)
+
+if(LeeYu && effects == "spfe"){
+ T <- T- 1
+ NT <- n*T
+}
+
+if(LeeYu && effects == "tpfe"){
+ n <- n-1
+ NT <- n*T
+}
+
+ if(LeeYu && effects == "sptpfe"){
+ n <- n-1
+ T <- T-1
+ NT <- n*T
+}
+ one <- T*(tr(WA %*% WA) + tr(t(WA) %*% WA))
two <- 1/as.numeric(s2) * t(betas) %*% xtWAWAxt %*% betas
- V <- one + two
- zero <- rbind(rep(0, length(betas)))
+ V <- one + two
+ zero <- rbind(rep(0, length(betas)))
col1 <- rbind(NT/(2 * (s2^2)), T*tr(WA)/s2, t(zero))
three <- (1/as.numeric(s2)) * xtWAxt %*% betas
col2 <- rbind(T*tr(WA)/s2, V, three )
@@ -141,16 +193,33 @@
-f_sarpanel_hess <- function (coefs, env)
+f_sarpanel_hess <- function (coefs, env, effects = effects, LeeYu = LeeYu)
{
T<-get("T", envir = env)
NT<-get("NT", envir = env)
+ n <- NT/T
+
+if(LeeYu && effects == "spfe"){
+ T <- T- 1
+ NT <- n*T
+}
+
+if(LeeYu && effects == "tpfe"){
+ n <- n-1
+ NT <- n*T
+}
+
+
+if(LeeYu && effects == "sptpfe"){
+ n <- n-1
+ T <- T-1
+ NT <- n*T
+}
+
lambda <- coefs[1]
beta <- coefs[-1]
-
- n <- NT/T
SSE <- sar_hess_sse_panel(lambda, beta, env)
s2 <- SSE /n
@@ -205,7 +274,7 @@
-sperrorlm<-function(env, zero.policy = zero.policy, interval = interval, con = con, llprof = llprof, tol.solve= tol.solve, Hess = Hess){
+sperrorlm <- function(env, zero.policy = zero.policy, interval = interval, con = con, llprof = llprof, tol.solve= tol.solve, Hess = Hess, LeeYu = LeeYu, effects = effects){
xt <- get("xt", envir = env)
yt <- get("yt", envir = env)
@@ -218,6 +287,7 @@
con<-get("con", envir = env)
NT<-get("NT", envir = env)
T<-get("T", envir = env)
+n <- NT/T
listw<-get("listw", envir = env)
inde<-get("inde", envir = env)
interval <- get("interval1", envir = env)
@@ -239,18 +309,62 @@
1)
r <- as.vector(residuals(lm.target))
p <- lm.target$rank
- s2 <- crossprod(r)/NT
+ SSE <- crossprod(residuals(lm.target))
+ s2 <- as.numeric(SSE)/NT
+
+if(LeeYu && effects == "spfe") s2 <- (T/(T-1)) * as.numeric(s2)
+if(LeeYu && effects == "tpfe") s2 <- (n/(n-1)) * as.numeric(s2)
+
rest.se <- (summary(lm.target)$coefficients[, 2]) * sqrt((NT - p)/NT)
betas <- coefficients(lm.target)
names(betas) <- colnames(xt)
coefs <- c(rho, betas)
+if(LeeYu && effects == "sptpfe"){
+
+ tr <- function(A) sum(diag(A))
+ W <- listw2dgCMatrix(listw, zero.policy = zero.policy)
+ A <- solve(sparseMatrix(i=1:(NT/T), j=1:(NT/T), x=1) - rho * W)
+ WA <- W %*% A
+ asyvar <- matrix(0, nrow = 2 + p, ncol = 2 + p)
+ asyvar[1, 1] <- NT/(2 * (s2^2))
+ asyvar[2, 1] <- asyvar[1, 2] <- T*tr(WA)/s2
+ asyvar[2, 2] <- T*(tr(WA %*% WA) + tr(t(WA) %*% WA))
+ asyvar[3:(p + 2), 3:(p + 2)] <- 1/as.numeric(s2) * (t(xt - rho *wxt) %*% (xt - rho * wxt))
+ asyv <- solve(asyvar, tol = con$tol.solve)
+ rownames(asyv) <- colnames(asyv) <- c("sigma","rho", colnames(xt))
+ s2.se <- sqrt(asyv[1, 1])
+ rho.se <- asyv[2, 2]
+ asyvar1 <- asyv[-c(1,2),-c(1,2)]
+ init <- c((T/(T+1)), rep(1,p+1))
+
+ a3 <- rep(0,p)
+ a2 <- 1/(1 - rho)
+ a1 <- 1/(2*s2)
+ a <- c(a1,a2,a3)
+ Bhat <- - (asyv/n) %*% a
+
+
+ coefs1 <- c(s2, rho, betas)
+ Theta2 <- init * coefs1 + Bhat
+ betas <- as.numeric(Theta2[(3:(2+p))])
+ names(betas) <- colnames(xt)
+ rho <-Theta2[2]
+ names(rho) <- "rho"
+ s2 <- Theta2[1]
+ coefs <- c(rho, betas)
+
+
+
+
+}
+
if(Hess){
- fd <- fdHess(coefs, sarpanelerror_hess, env)
+ fd <- fdHess(coefs, sarpanelerror_hess, env, LeeYu = LeeYu, effects = effects)
mat <- fd$Hessian
- fdHess<- solve(-(mat), tol.solve = tol.solve)
+ fdHess<- solve(-(mat), tol.solve = tol.solve)
rownames(fdHess) <- colnames(fdHess) <- c("rho",colnames(xt))
rho.se <- fdHess[1, 1]
@@ -263,6 +377,27 @@
W <- listw2dgCMatrix(listw, zero.policy = zero.policy)
A <- solve(sparseMatrix(i=1:(NT/T), j=1:(NT/T), x=1) - rho * W)
WA <- W %*% A
+
+if(LeeYu && effects == "spfe"){
+
+ T <- T- 1
+ NT <- n*T
+
+}
+
+if(LeeYu && effects == "tpfe"){
+
+ n <- n-1
+ NT <- n*T
+
+}
+
+ if(LeeYu && effects == "sptpfe"){
+ n <- n-1
+ T <- T-1
+ NT <- n*T
+}
+
asyvar <- matrix(0, nrow = 2 + p, ncol = 2 + p)
asyvar[1, 1] <- NT/(2 * (s2^2))
asyvar[2, 1] <- asyvar[1, 2] <- T*tr(WA)/s2
@@ -282,7 +417,7 @@
-sarpanelerror_hess<-function (coef, env=env)
+sarpanelerror_hess<-function (coef, env=env, LeeYu = LeeYu, effects = effects)
{
yt<- get("yt", envir = env)
xt<- get("xt", envir = env)
@@ -295,16 +430,29 @@
NT<- get("NT", envir = env)
inde<- get("inde", envir = env)
T<- get("T", envir = env)
+ n <- NT/T
-# coefsl <- c(s2, lambda, betas)
+if(LeeYu && effects == "spfe"){
+ T <- T- 1
+ NT <- n*T
+}
- # s2 <- coef[1]
+if(LeeYu && effects == "tpfe"){
+ n <- n-1
+ NT <- n*T
+}
+
+ if(LeeYu && effects == "sptpfe"){
+ n <- n-1
+ T <- T-1
+ NT <- n*T
+}
+
rho <- coef[1]
bb <- coef[-1]
yco <- yt - rho * wyt
xco <- xt - rho * wxt
- # bb<- solve(crossprod(xco),crossprod(xco, yco) )
ehat<- yco - xco %*% bb
SSE <- crossprod(ehat)
@@ -321,18 +469,13 @@
###SARAR MODEL
-sacsarpanel<-function (coefs, env, LeeYu = LeeYu){
+sacsarpanel<-function (coefs, env){
lambda <- coefs[1]
rho <- coefs[2]
T<-get("T", envir = env)
n <- get("n", envir = env)
-if(LeeYu){
- T <- T-1
- NT <- n*T
- }
-
SSE <- sacsarpanel_sse(coefs, env)
s2 <- SSE/n
ldet1 <- do_ldet(lambda, env, which = 1)
@@ -340,10 +483,9 @@
ret <- (T * ldet1 + T * ldet2 - (((n*T)/2) * (log(2 * pi)+1)) - (n*T/2) * log(s2))
# - (1/(2 * (s2))) * SSE)
-
- # if(get("verbose", envir = env)) cat("lambda:", lambda, " rho:", rho, " function:",
- # ret, " Jacobian1:", ldet1, " Jacobian2:", ldet2,
- # " SSE:", SSE, "\n")
+if(get("verbose", envir = env)) cat("lambda:", lambda, " rho:", rho, " function:",
+ ret, " Jacobian1:", ldet1, " Jacobian2:", ldet2,
+ " SSE:", SSE, "\n")
-ret
}
@@ -363,7 +505,7 @@
}
-spsararlm<-function(env, zero.policy = zero.policy, con = con, llprof = llprof, tol.solve= tol.solve, Hess = Hess, method = method, LeeYu = LeeYu){
+spsararlm<-function(env, zero.policy = zero.policy, con = con, llprof = llprof, tol.solve= tol.solve, Hess = Hess, method = method, LeeYu = LeeYu, effects = effects){
xt <- get("xt", envir = env)
yt <- get("yt", envir = env)
@@ -379,11 +521,6 @@
T<-get("T", envir = env)
n<-get("n", envir = env)
-if(LeeYu){
- T <- T-1
- NT <- n*T
- }
-
listw<-get("listw", envir = env)
listw2<-get("listw2", envir = env)
inde<-get("inde", envir = env)
@@ -430,42 +567,42 @@
if (is.null(pars)) {
mxs <- apply(mpars, 1, function(pp) -nlminb(pp, sacsarpanel,
env = env, control = con$opt_control, lower = lower,
- upper = upper, LeeYu = LeeYu)$objective)
+ upper = upper)$objective)
pars <- mpars[which.max(mxs), ]
optres <- nlminb(pars, sacsarpanel, env = env, control = con$opt_control,
- lower = lower, upper = upper, LeeYu = LeeYu)
+ lower = lower, upper = upper)
}
else {
optres <- nlminb(pars, sacsarpanel, env = env, control = con$opt_control,
- lower = lower, upper = upper, LeeYu = LeeYu)
+ lower = lower, upper = upper)
}
}
else if (con$opt_method == "L-BFGS-B") {
if (is.null(pars)) {
mxs <- apply(mpars, 1, function(pp) -optim(pars,
sacsarpanel, env = env, method = "L-BFGS-B", control = con$opt_control,
- lower = lower, upper = upper, LeeYu = LeeYu)$objective)
+ lower = lower, upper = upper)$objective)
pars <- mpars[which.max(mxs), ]
optres <- optim(pars, sacsarpanel, env = env, method = "L-BFGS-B",
- control = con$opt_control, lower = lower, upper = upper, LeeYu = LeeYu)
+ control = con$opt_control, lower = lower, upper = upper)
}
else {
optres <- optim(pars, sacsarpanel, env = env, method = "L-BFGS-B",
- control = con$opt_control, lower = lower, upper = upper, LeeYu = LeeYu)
+ control = con$opt_control, lower = lower, upper = upper)
}
}
else {
if (is.null(pars)) {
mxs <- apply(mpars, 1, function(pp) -optim(pars,
sacsarpanel, env = env, method = con$opt_method,
- control = con$opt_control, LeeYu = LeeYu)$objective)
+ control = con$opt_control)$objective)
pars <- mpars[which.max(mxs), ]
optres <- optim(pars, sacsarpanel, env = env, method = con$opt_method,
- control = con$opt_control, LeeYu = LeeYu)
+ control = con$opt_control)
}
else {
optres <- optim(pars, sacsarpanel, env = env, method = con$opt_method,
- control = con$opt_control, LeeYu = LeeYu)
+ control = con$opt_control)
}
}
@@ -479,9 +616,6 @@
lambda <- optres$par[1]
names(lambda) <- "lambda"
- # timings[["coefs"]] <- proc.time() - .ptime_start
- # .ptime_start <- proc.time()
- # assign("first_time", TRUE, envir = env)
lm.target <- lm(I(yt - lambda * wyt - rho * w2yt + rho * lambda *
@@ -490,17 +624,101 @@
r <- as.vector(residuals(lm.target))
fit <- as.vector(yt - r)
p <- lm.target$rank
- SSE <- deviance(lm.target)
- s2 <- SSE/NT
+ SSE <- crossprod(residuals(lm.target))
+ s2 <- as.numeric(SSE)/NT
+
+if(LeeYu && effects == "spfe") s2 <- (T/(T-1)) * as.numeric(s2)
+if(LeeYu && effects == "tpfe") s2 <- (n/(n-1)) * as.numeric(s2)
+
betas <- coefficients(lm.target)
names(betas) <- colnames(xt)
- # coefs <- c(rho, lambda, betas)
coefs <- c(lambda, rho, betas)
+if(LeeYu && effects == "sptpfe"){
+
+ tr <- function(A) sum(diag(A))
+ W1 <- listw2dgCMatrix(listw, zero.policy = zero.policy)
+ W2 <- listw2dgCMatrix(listw2, zero.policy = zero.policy)
+ Sl <- sparseMatrix(i=1:(NT/T), j=1:(NT/T), x=1) - lambda * W1
+ Rr <- sparseMatrix(i=1:(NT/T), j=1:(NT/T), x=1) - rho * W2
+ Slinv <- solve(Sl)
+ Rrinv <- solve(Rr)
+ Gmat <- W1 %*% Slinv
+ Hmat <- W2 %*% Rrinv
+ It <- sparseMatrix(i=1:T, j=1:T, x=1)
+ Jn <- Diagonal(n) - (1/n) * outer(rep(1,n),rep(1,n))
+
+
+# Equation 53 Lee and Yu
+ Wdot <- Rr %*% W1 %*% Rrinv
+ Gdot <- Rr %*% Gmat %*% Rrinv
+ GS <- t(Gdot) + Gdot
+ HS <- t(Hmat) + Hmat
+ Rrbig <- kronecker(It,Rr)
+ RriB <- kronecker(It,Rrinv)
+ GdotB<- kronecker(It,Gdot)
+ WdotB<- kronecker(It,Wdot)
+ JnB <- kronecker(It,Jn)
+ Xdot <- Rrbig %*% xt
+ JXdot <- JnB %*% Xdot
+ GdXdb <- GdotB %*% Xdot %*% betas
+ JGdXdb <- JnB %*% GdotB %*% Xdot %*% betas
+
+ fp <- (1/s2) * crossprod(GdXdb, JGdXdb)
+ lala <- fp + (T * tr(GS %*% Jn %*% Gdot))
+ laro <- T * tr(HS %*% Jn %*% Gdot)
+ lasi <- (1/s2) * T * tr(Gdot)
+ roro <- T * tr(HS %*% Hmat)
+ rosi <- (1/s2) * T * tr(Hmat)
+ sisi <- NT/(2*s2*s2)
+ bebe <- (1/s2) * crossprod(Xdot, JXdot)
+ bela <- (1/s2) * crossprod(GdXdb, JXdot)
+
+ asyvar <- matrix(0, nrow = 3 + p, ncol = 3 + p)
+ asyvar[1:p, 1:p] <- as.matrix(bebe)
+ asyvar[p+1, 1:p] <- asyvar[1:p, p+1] <- as.numeric(bela)
+ asyvar[p+2, 1:p] <- asyvar[1:p, p+2] <- 0
+ asyvar[p+3, 1:p] <- asyvar[1:p, p+3] <- 0
+ asyvar[p+2, p+1] <- asyvar[p+1, p+2] <- as.numeric(laro)
+ asyvar[p+3, p+1] <- asyvar[p+1, p+3] <- as.numeric(lasi)
+ asyvar[p+3, p+2] <- asyvar[p+2, p+3] <- as.numeric(rosi)
+ asyvar[1+p, 1+p] <- as.matrix(lala)
+ asyvar[2+p, 2+p] <- as.matrix(roro)
+ asyvar[3+p, 3+p] <- as.matrix(sisi)
+ asyv <- solve(asyvar, tol = con$tol.solve)
+
+ a1 <- rep(0,p)
+ a2 <- as.numeric((1/n) * rep(1,n) %*% Gdot %*% rep(1,n))
+ a3 <- as.numeric((1/n) * rep(1,n) %*% Hmat %*% rep(1,n))
+ a4 <- as.numeric(1/(2*s2))
+ a <- c(a1,a2,a3,a4)
+ Bhat <- - asyv %*% a
+ At <- matrix(0, nrow = 3 + p, ncol = 3 + p)
+ At[(1:(p+2)), (1:(p+2))]<- diag((p+2))
+ At[(1:(p+2)), 3+p] <- rep(0,p+2)
+ At[3+p,(1:(p+2))] <- rep(0,p+2)
+ At[3+p, 3+p] <- T/(T-1)
+ coefs1 <- c(betas, lambda, rho, s2)
+ Theta1 <- coefs1 - (Bhat/n)
+ Theta2 <- At %*% Theta1
+ betas <- Theta2[1:p]
+ names(betas) <- colnames(xt)
+ lambda <-Theta2[p+1]
+ rho <- Theta2[p+2]
+ names(rho) <- "rho"
+ names(lambda) <- "lambda"
+ s2 <- Theta2[p+3]
+ coefs <- c(lambda, rho, betas)
+
+}
+
+
###Add the vc matrix exact
-if(Hess){
- fd <- fdHess(coefs, f_sacpanel_hess, env, LeeYu = LeeYu)
- #
+if(Hess){
+
+# if(LeeYu && effects == "sptpfe") stop("Numerical Hessian should not be calculated when 'LeeYu = TRUE' and effects are 'twoways' ")
+
+ fd <- fdHess(coefs, f_sacpanel_hess, env, LeeYu = LeeYu, effects = effects)
mat <- fd$Hessian
fdHess<- solve(-(mat), tol.solve = tol.solve)
rownames(fdHess) <- colnames(fdHess) <- c("lambda", "rho",colnames(xt))
@@ -512,7 +730,8 @@
}
else{
-
+
+
tr <- function(A) sum(diag(A))
W1 <- listw2dgCMatrix(listw, zero.policy = zero.policy)
W2 <- listw2dgCMatrix(listw2, zero.policy = zero.policy)
@@ -522,22 +741,8 @@
Rrinv <- solve(Rr)
Gmat <- W1 %*% Slinv
Hmat <- W2 %*% Rrinv
- # Gmat2 <- Gmat %*% Gmat
- # Hmat2 <- Hmat %*% Hmat
+ It <- sparseMatrix(i=1:T, j=1:T, x=1)
-if(LeeYu) dim <- T+1
-else dim <- T
- It <- sparseMatrix(i=1:dim, j=1:dim, x=1)
-
-
- # Slinv <- kronecker(It,Slinv)
- # Rrinv <- kronecker(It,Rrinv)
- # Gmat <- kronecker(It,Gmat)
- # Hmat <- kronecker(It,Hmat)
- # W1 <- kronecker(It,W1)
- # W2 <- kronecker(It,W2)
- # Rr<- kronecker(It,Rr)
- # Sl<- kronecker(It,Sl)
# Equation 39 Lee and Yu
Wdot <- Rr %*% W1 %*% Rrinv
@@ -550,9 +755,23 @@
WdotB<- kronecker(It,Wdot)
Xdot <- Rrbig %*% xt
GdXdb <- GdotB %*% Xdot %*% betas
+
+if(LeeYu && effects == "spfe"){
+ T <- T- 1
+ NT <- n*T
+}
+if(LeeYu && effects == "tpfe"){
+ n <- n-1
+ NT <- n*T
+}
+if(LeeYu && effects == "sptpfe"){
+ n <- n-1
+ T <- T-1
+ NT <- n*T
+}
fp <- (1/s2) *crossprod(GdXdb)
lala <- fp + (T * tr(GS %*% Gdot))
laro <- T * tr(HS %*% Gdot)
@@ -563,36 +782,7 @@
bebe <- (1/s2) * crossprod(Xdot)
bela <- (1/s2) * crossprod(GdXdb, Xdot)
- # Rrxt <- Rr %*% xt
- # bebe <- (1/s2) * crossprod(Rrxt)
- # RrWyt <- Rr %*% W1 %*% yt
- # bela <- (1/s2) * crossprod(RrWyt, Rrxt)
- # Vxsi <- Rr %*% (Sl %*% yt - xt %*% betas)
- # HVxsi <- Hmat %*% Vxsi
- # bero1 <- (1/s2) * crossprod(HVxsi, Rrxt)
- # Mxt <- W2 %*% xt
- # bero2 <- (1/s2) * crossprod(Vxsi, Mxt)
- # bero <- bero1 + bero2
- # besi<- (1/s2) * (1/s2) * crossprod(Vxsi, Rrxt)
- # lala1 <- (1/s2) *crossprod(RrWyt)
- # lala2 <- T* tr(Gmat2)
- # lala <- lala1 + lala2
- # laro1 <- (1/s2) *crossprod(RrWyt, HVxsi)
- # MWyt <- W2 %*% W1 %*% yt
- # laro2 <- (1/s2) *crossprod(MWyt, Vxsi)
- # laro <- laro1 + laro2
- # lasi <- (1/s2) * (1/s2) * crossprod(RrWyt, Vxsi)
- # roro1 <- (1/s2) * crossprod(HVxsi)
- # roro2 <- T * tr(Hmat2)
- # roro <- roro1 + roro2
- # rosi <- (1/s2) * (1/s2) * crossprod(HVxsi, Vxsi)
- # sisi1 <- - (n*T)/(2*s2*s2)
- # sisi2 <- (1/(s2)^3) * crossprod(Vxsi)
- # sisi <- sisi1 + sisi2
-
asyvar <- matrix(0, nrow = 3 + p, ncol = 3 + p)
- # print(dim(bela))
- # print(p)
asyvar[1:p, 1:p] <- as.matrix(bebe)
asyvar[p+1, 1:p] <- asyvar[1:p, p+1] <- as.numeric(bela)
asyvar[p+2, 1:p] <- asyvar[1:p, p+2] <- 0
@@ -603,37 +793,41 @@
asyvar[1+p, 1+p] <- as.matrix(lala)
asyvar[2+p, 2+p] <- as.matrix(roro)
asyvar[3+p, 3+p] <- as.matrix(sisi)
-# print(asyvar)
-
- # asyvar <- NT * asyvar
asyv <- solve(asyvar, tol = con$tol.solve)
- # print(sqrt(diag(asyv)))
rownames(asyv) <- colnames(asyv) <- c(colnames(xt), "lambda", "rho", "sigma")
s2.se <- asyv[3+p, 3+p]
rho.se <- asyv[2+p, 2+p]
lambda.se <- asyv[1+p, 1+p]
- # print(rho.se)
- # print(lambda.se)
rest.se <- sqrt(diag(asyv))[-((p+1):(p+3))]
asyvar1 <- asyv[-((p+1):(p+3)),-((p+1):(p+3))]
-# stop("Asymptotic VC matrix not yet implemented for model SARAR")
}
return<-list(coeff = betas, lambda = lambda, rho = rho, s2 = s2, asyvar1 = asyvar1, lambda.se = lambda.se, rho.se = rho.se, s2.se = s2.se)
}
-f_sacpanel_hess <- function (coefs, env, LeeYu = LeeYu)
+f_sacpanel_hess <- function (coefs, env, LeeYu = LeeYu, effects = effects)
{
T<-get("T", envir = env)
NT<-get("NT", envir = env)
n<-get("n", envir = env)
-if(LeeYu){
+if(LeeYu && effects == "spfe"){
+ T <- T- 1
+ NT <- n*T
+}
+
+if(LeeYu && effects == "tpfe"){
+ n <- n-1
+ NT <- n*T
+}
+
+if(LeeYu && effects == "sptpfe"){
+ n <- n-1
T <- T-1
NT <- n*T
- }
+}
lambda <- coefs[1]
rho <- coefs[2]
Modified: pkg/R/spfeml.R
===================================================================
--- pkg/R/spfeml.R 2013-11-06 23:12:32 UTC (rev 179)
+++ pkg/R/spfeml.R 2013-11-23 22:56:34 UTC (rev 180)
@@ -1,11 +1,13 @@
-spfeml<-function(formula, data=list(), index=NULL, listw, listw2 = NULL, na.action, model = c("lag","error", "sarar"),effects = c('spfe','tpfe','sptpfe'), method="eigen", quiet = TRUE, zero.policy = NULL, interval1 = NULL, interval2 = NULL, trs1 = NULL, trs2 = NULL, tol.solve = 1e-10, control = list(), legacy = FALSE, llprof = NULL, cl = NULL, Hess = TRUE, LeeYu = FALSE, ...){
+spfeml<-function(formula, data=list(), index=NULL, listw, listw2 = NULL, na.action, model = c("lag","error", "sarar"), effects = c('spfe','tpfe','sptpfe'), method="eigen", quiet = TRUE, zero.policy = NULL, interval1 = NULL, interval2 = NULL, trs1 = NULL, trs2 = NULL, tol.solve = 1e-10, control = list(), legacy = FALSE, llprof = NULL, cl = NULL, Hess = TRUE, LeeYu = FALSE, ...){
# timings <- list()
# .ptime_start <- proc.time()
model<-match.arg(model)
+effects <- match.arg(effects)
+
if (model == "sarar") con <- list(LAPACK = FALSE, Imult = 2L, cheb_q = 5L, MC_p = 16L, MC_m=30L, super=FALSE, opt_method = "nlminb", opt_control = list(), pars = NULL, npars = 4L, pre_eig1 = NULL, pre_eig2 = NULL)
else con <- list(tol.opt = .Machine$double.eps^0.5, Imult = 2, cheb_q = 5, MC_p = 16, MC_m = 30, super = NULL, spamPivot = "MMD", in_coef = 0.1, type = "MC", correct = TRUE, trunc = TRUE, SE_method = "LU", nrho = 200, interpn = 2000, SElndet = NULL, LU_order = FALSE, pre_eig = NULL)
@@ -315,66 +317,15 @@
assign("con", con, envir=env)
-# timings[["set_up"]] <- proc.time() - .ptime_start
-# .ptime_start <- proc.time()
-
-
-#Lee and Yu transformation
-
-# if(LeeYu){
-# T <- T-1
-# assign("T",T, envir=env)
-# NT <- n*T
-# assign("NT",T, envir=env)
-# # stop("Lee and Yu correction not yet implemented")
-# # if (effects=="spfe" | effects=="sptpfe"){
-# # IT <- Diagonal(T)
-# # IN <- Diagonal(n)
-# # JT <- matrix(1,T,T)
-# # Jbar <- 1/T * JT
-# # Qmat <-IT - Jbar
-# # vec <- eigen(Qmat)
-# # Fmat <- vec$vectors[,vec$values==1L]
-# # Ftm <- kronecker(t(Fmat), IN)
-# # }
-
-
-# # if (effects=="tpfe" | effects=="sptpfe"){
-# # IT <- Diagonal(T)
-# # IN <- Diagonal(n)
-# # JT <- matrix(1,T,T)
-# # Jbar <- 1/T * JT
-# # Qmat <-IT - Jbar
-# # vec <- eigen(Qmat)
-# # Fmat <- matrix(vec$vectors[,vec$values==1L], T, T-1)
-# # Ftm <- kronecker(t(Fmat), IN)
-# # iotan <- matrix(1,n,1)
-# # Jnbar <-1/n * iotan %*% t(iotan)
-# # Qmat1 <- IN - Jnbar
-# # vec1 <- eigen(Qmat1)
-# # Fmat1 <- matrix(vec1$vectors[,vec1$values==1L], n, n-1)
-# # FFmat<- kronecker(t(Fmat), t(Fmat1))
-# # }
-
-
-
-# }
-
-
-
-
if (!quiet)
cat(paste("\nSpatial fixed effects model\n", "Jacobian calculated using "))
if(model == "lag"){
interval1 <- spdep:::jacobianSetup(method, env, con, pre_eig = con$pre_eig, trs = trs1, interval = interval1)
assign("interval1", interval1, envir = env)
- # nm <- paste(method, "set_up", sep = "_")
- # timings[[nm]] <- proc.time() - .ptime_start
- # .ptime_start <- proc.time()
- RES<- splaglm(env = env, zero.policy = zero.policy, interval = interval1, Hess = Hess)
+ RES<- splaglm(env = env, zero.policy = zero.policy, interval = interval1, con = con, llprof = llprof, tol.solve= tol.solve, Hess = Hess, method = method, LeeYu = LeeYu, effects = effects)
res.eff<-felag(env = env, beta = RES$coeff, sige = RES$s2, effects = effects, method = method, lambda = RES$lambda, legacy = legacy, zero.policy = zero.policy)
@@ -390,7 +341,7 @@
# timings[[nm]] <- proc.time() - .ptime_start
# .ptime_start <- proc.time()
- RES<- spsararlm(env = env, zero.policy = zero.policy, con = con, llprof = llprof, tol.solve = tol.solve, Hess = Hess, LeeYu = LeeYu)
+ RES<- spsararlm(env = env, zero.policy = zero.policy, con = con, llprof = llprof, tol.solve = tol.solve, Hess = Hess, LeeYu = LeeYu, effects = effects)
res.eff<-felag(env = env, beta=RES$coeff, sige=RES$s2, effects = effects ,method = method, lambda = RES$lambda, legacy = legacy, zero.policy = zero.policy)
@@ -407,7 +358,7 @@
# timings[[nm]] <- proc.time() - .ptime_start
# .ptime_start <- proc.time()
- RES<- sperrorlm(env = env, zero.policy = zero.policy, interval = interval1, Hess = Hess)
+ RES<- sperrorlm(env = env, zero.policy = zero.policy, interval = interval1, Hess = Hess, LeeYu = LeeYu, effects = effects)
res.eff<-feerror(env = env, beta=RES$coeff, sige=RES$s2, effects = effects ,method =method, rho=RES$rho, legacy = legacy)
}
Modified: pkg/R/spgm.R
===================================================================
--- pkg/R/spgm.R 2013-11-06 23:12:32 UTC (rev 179)
+++ pkg/R/spgm.R 2013-11-23 22:56:34 UTC (rev 180)
@@ -127,7 +127,7 @@
if (instr) H <- Hinst
else H <- cbind(X, Zinst)
Z <- cbind(yend, X)
- Znames<- colnames(Z)
+ Znames<- colnames(Z)
yendp<-matrix(,nrow(yend), ncol(yend))
for(i in 1:ncol(yend)) yendp[,i] <- fitted.values(lm(yend[,i] ~ H-1 ))
Zp <- cbind(yendp,X)
@@ -135,7 +135,7 @@
biv <- coefficients(model.fit)
readout<- which(is.na(biv))
-if(any(is.na(biv))) yp <- Z[,-which(is.na(biv))] %*% biv[-which(is.na(biv))]
+if(any(is.na(biv))) yp <- as.matrix(Z[,-which(is.na(biv))]) %*% biv[-which(is.na(biv))]
else yp <- Z %*% biv
@@ -145,7 +145,7 @@
df <- model.fit$df.residual
s2 <- sse/df
-if(any(is.na(biv))) Zp<-Zp[,-which(is.na(biv))]
+if(any(is.na(biv))) Zp<-as.matrix(Zp)[,-which(is.na(biv))]
ZpZi<-solve(crossprod(Zp))
@@ -681,6 +681,7 @@
y <- model.extract(mf, "response")
x <- model.matrix(mt, mf)
+ namesx <- colnames(x)
N<-length(unique(ind))
k<-dim(x)[[2]]
@@ -715,7 +716,7 @@
Xbetween <- as.matrix(Xbetween)
Xwithin <- as.matrix(Xwithin)
del<-which(diag(var(as.matrix(Xwithin)))==0)
-if (colnames(x)[1] == "(Intercept)") Xbetween <- Xbetween[,-1]
+if (namesx[1] == "(Intercept)") Xbetween <- Xbetween[,-1]
delb<-which(diag(var(as.matrix(Xbetween)))==0)
if(length(delb)==0) Xbetween<-Xbetween
else Xbetween<-Xbetween[,-delb]
@@ -757,7 +758,7 @@
finrho<-estim1$par[1]
finsigmaV<-estim1$par[2]
-#print(c(finrho,finsigmaV))
+# print(c(finrho,finsigmaV))
wy <- as.matrix(Ws2 %*% y)
yt <- y-finrho*wy
@@ -773,7 +774,7 @@
wyf<-panel.transformations(wyt,indic, type= "within")
xf<-xf[,-del]
xf<-as.matrix(xf)
- colnames(xf)<-colnames(x)[-del]
+ colnames(xf)<- namesx[-del]
# wxf <- as.matrix(Ws %*% xf)
Modified: pkg/R/spml.R
===================================================================
--- pkg/R/spml.R 2013-11-06 23:12:32 UTC (rev 179)
+++ pkg/R/spml.R 2013-11-23 22:56:34 UTC (rev 180)
@@ -11,21 +11,36 @@
## check class(listw)
checklw <- function(x) {
+# if(model == "within"){
+
+ # if("matrix" %in% class(x)) x <- Matrix(x)
+ # if("listw" %in% class(x)) x <- listw2dgCMatrix(x)
+ # if("Matrix" %in% class(x)) x <- x
+
+ # else stop("'listw' has to be 'listw', 'matrix', or 'Matrix' when model is within")
+
+ # }
+# else{
+
if(!("listw" %in% class(x))) {
if("matrix" %in% class(x)) {
require(spdep)
x <- listw2mat(x)
- } else {
+ }
+ else {
stop("'listw' has to be either a 'listw' or a 'matrix' object")
}}
+ # }
return(x)
}
+
checklw(listw)
checklw(listw2)
## dimensions check is moved downstream
switch(match.arg(model), within={
+
if(lag) {
model <- switch(match.arg(spatial.error), b="sarar",
kkp="sarar", none="lag")
More information about the Splm-commits
mailing list