[Splm-commits] r161 - pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Apr 5 01:34:29 CEST 2013
Author: gpiras
Date: 2013-04-05 01:34:29 +0200 (Fri, 05 Apr 2013)
New Revision: 161
Modified:
pkg/R/likelihoodsFE.R
Log:
updated spfeml
Modified: pkg/R/likelihoodsFE.R
===================================================================
--- pkg/R/likelihoodsFE.R 2013-04-02 21:02:16 UTC (rev 160)
+++ pkg/R/likelihoodsFE.R 2013-04-04 23:34:29 UTC (rev 161)
@@ -327,7 +327,8 @@
rho <- coefs[2]
T<-get("T", envir = env)
n <- get("n", envir = env)
- if(LeeYu){
+
+if(LeeYu){
T <- T-1
NT <- n*T
}
@@ -505,74 +506,103 @@
rho.se <- fdHess[2,2]
lambda.se <- fdHess[1,1]
asyvar1 <- vcov(lm.target)
+ s2.se <- NULL
}
else{
- tr <- function(A) sum(diag(A))
- W1 <- listw2dgCMatrix(listw, zero.policy = zero.policy)
- W2 <- listw2dgCMatrix(listw2, zero.policy = zero.policy)
- Sl <- solve(sparseMatrix(i=1:(NT/T), j=1:(NT/T), x=1) - lambda * W1)
- Rr <- solve(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
- Rrxt <- Rr %*% xt
- Gmat2 <- Gmat %*% Gmat
- Hmat2 <- Hmat %*% Hmat
- bebe <- (1/s2) * crossprod(Rrxt)
- RrWyt <- Rr %*% W1 %*% yt
- bela <- (1/s2) * crossprod(RrWyt, Rrxt)
- Vxsi <- Rr %*% (Sl %*% yt - xt %*% beta)
- 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 <- ll1 + ll2
- 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) * (1/s2) * (1/s2) * crossprod(Vxsi)
- sisi <- sisi1 + sisi2
+ # tr <- function(A) sum(diag(A))
+ # W1 <- listw2dgCMatrix(listw, zero.policy = zero.policy)
+ # W2 <- listw2dgCMatrix(listw2, zero.policy = zero.policy)
+ # Sl <- solve(sparseMatrix(i=1:(NT/T), j=1:(NT/T), x=1) - lambda * W1)
+ # Rr <- solve(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
+ # Gmat2 <- Gmat %*% Gmat
+ # Hmat2 <- Hmat %*% Hmat
+ # It <- sparseMatrix(i=1:(T+1), j=1:(T+1), x=1)
- asyvar <- matrix(0, nrow = 3 + p, ncol = 3 + p)
- asyvar[1:p, 1:p] <- bebe
- asyvar[p+1, 1] <- asyvar[1, p+1] <- bela
- asyvar[p+2, 1] <- asyvar[1, p+2] <- bero
- asyvar[p+3, 1] <- asyvar[1, p+3] <- besi
- asyvar[p+2, p+1] <- asyvar[p+1, p+2] <- laro
- asyvar[p+3, p+1] <- asyvar[p+1, p+3] <- lasi
- asyvar[p+3, p+2] <- asyvar[p+2, p+3] <- rosi
- asyvar[1+p, 1+p] <- lala
- asyvar[2+p, 2+p] <- roro
- asyvar[3+p, 3+p] <- sisi
+
+ # 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)
- asyv <- solve(asyvar, tol = con$tol.solve)
- rownames(asyv) <- colnames(asyv) <- c(colnames(xt), "lambda", "rho", "sigma")
- s2.se <- sqrt(asyv[3+p, 3+p])
- rho.se <- asyv[2+p, 2+p]
- lambda.se <- asyv[2+p, 2+p]
- asyvar1 <- asyv[-((p+1):(p+3)),-((p+1):(p+3))]
+# Equation 43 Lee and Yu
+ # Wdot <- Rr %*% W1 %*% Rrinv
+ # Gdot <- Rr %*% Gmat %*% Rrinv
+ # Rrbig <- kronecker(It,Rr)
+ # GdotB<- kronecker(It,Gdot)
+ # WdotB<- kronecker(It,Wdot)
+
+
+ # Xdot <- Rrbig %*% xt
+
+
+
+ # 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] <- as.numeric(bero)
+ # asyvar[p+3, 1:p] <- asyvar[1:p, p+3] <- as.numeric(besi)
+ # 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)
+# print(asyvar)
+ # asyv <- solve(asyvar, tol = con$tol.solve)
+# # print(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")
+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)
+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)
More information about the Splm-commits
mailing list