[Gmm-commits] r13 - pkg/gmm/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Dec 21 02:10:45 CET 2009
Author: chaussep
Date: 2009-12-21 02:10:45 +0100 (Mon, 21 Dec 2009)
New Revision: 13
Modified:
pkg/gmm/R/Methods.gel.R
pkg/gmm/R/charStable.R
pkg/gmm/R/gel.R
pkg/gmm/R/specTest.R
Log:
reorganizing the codes
Modified: pkg/gmm/R/Methods.gel.R
===================================================================
--- pkg/gmm/R/Methods.gel.R 2009-12-18 01:23:31 UTC (rev 12)
+++ pkg/gmm/R/Methods.gel.R 2009-12-21 01:10:45 UTC (rev 13)
@@ -11,7 +11,7 @@
# A copy of the GNU General Public License is available at
# http://www.r-project.org/Licenses/
-confint.gel <- function(object, parm, level=0.95, lambda=FALSE, ...)
+confint.gel <- function(object, parm, level = 0.95, lambda = FALSE, ...)
{
z <- object
n <- nrow(z$gt)
@@ -23,26 +23,26 @@
se_parl <- sqrt(diag(z$vcov_lambda))
lamb <- z$lambda
- zs <- qnorm((1-level)/2,lower.tail=FALSE)
+ zs <- qnorm((1 - level)/2, lower.tail=FALSE)
ch <- zs*se_par
if(!lambda)
{
- ans <- cbind(par-ch,par+ch)
- dimnames(ans) <- list(names(par),c((1-level)/2,0.5+level/2))
+ ans <- cbind(par-ch, par+ch)
+ dimnames(ans) <- list(names(par), c((1 - level)/2, 0.5+level/2))
}
if(lambda)
{
chl <- zs*se_parl
- ans <- cbind(lamb-chl,lamb+chl)
- dimnames(ans) <- list(names(lamb),c((1-level)/2,0.5+level/2))
+ ans <- cbind(lamb - chl, lamb + chl)
+ dimnames(ans) <- list(names(lamb), c((1 - level)/2, 0.5 + level/2))
}
if(!missing(parm))
ans <- ans[parm,]
ans
}
-coef.gel <- function(object, lambda=FALSE, ...)
+coef.gel <- function(object, lambda = FALSE, ...)
{
if(!lambda)
object$coefficients
@@ -50,7 +50,7 @@
object$lambda
}
-vcov.gel <- function(object, lambda=FALSE, ...)
+vcov.gel <- function(object, lambda = FALSE, ...)
{
if(!lambda)
object$vcov_par
@@ -58,15 +58,15 @@
object$vcov_lambda
}
-print.gel <- function(x, digits=5, ...)
+print.gel <- function(x, digits = 5, ...)
{
- cat("Type de GEL: ", x$type,"\n\n")
+ cat("Type de GEL: ", x$type, "\n\n")
cat("Coefficients:\n")
- print.default(format(coef(x), digits=digits),
+ print.default(format(coef(x), digits = digits),
print.gap = 2, quote = FALSE)
cat("\n")
cat("Lambdas:\n")
- print.default(format(coef(x,lambda=TRUE), digits=digits),
+ print.default(format(coef(x, lambda = TRUE), digits = digits),
print.gap = 2, quote = FALSE)
invisible(x)
}
@@ -74,23 +74,23 @@
print.summary.gel <- function(x, digits = 5, ...)
{
cat("\nCall:\n")
- cat(paste(deparse(x$call), sep="\n", collapse = "\n"), "\n\n", sep="")
- cat("\nType of GEL: ", x$type,"\n\n")
- cat("Kernel: ", x$kernel,"\n\n")
+ cat(paste(deparse(x$call), sep = "\n", collapse = "\n"), "\n\n", sep = "")
+ cat("\nType of GEL: ", x$type, "\n\n")
+ cat("Kernel: ", x$kernel, "\n\n")
cat("Coefficients:\n")
- print.default(format(x$coefficients, digits=digits),
+ print.default(format(x$coefficients, digits = digits),
print.gap = 2, quote = FALSE)
cat("\nLambdas:\n")
print.default(format(x$lambda, digits=digits),
print.gap = 2, quote = FALSE)
- cat("\n",x$stest$ntest,"\n")
+ cat("\n", x$stest$ntest, "\n")
print.default(format(x$stest$test, digits=digits),
print.gap = 2, quote = FALSE)
- cat("\nConvergence code for the coefficients: ",x$conv_par,"\n")
- cat("\nConvergence code for the lambdas: ",x$conv_lambda,"\n")
+ cat("\nConvergence code for the coefficients: ", x$conv_par, "\n")
+ cat("\nConvergence code for the lambdas: ", x$conv_lambda, "\n")
invisible(x)
}
@@ -107,11 +107,11 @@
lamb <- z$lambda
tvall <- lamb/se_parl
- ans <- list(type=z$type,call=z$call)
+ ans <- list(type = z$type, call = z$call)
names(ans$type) <-"Type of GEL"
- ans$coefficients <- round(cbind(par,se_par, tval, 2 * pnorm(abs(tval), lower.tail = FALSE)),5)
- ans$lambda <- round(cbind(lamb,se_parl, tvall, 2 * pnorm(abs(tvall), lower.tail = FALSE)),5)
+ ans$coefficients <- round(cbind(par, se_par, tval, 2 * pnorm(abs(tval), lower.tail = FALSE)), 5)
+ ans$lambda <- round(cbind(lamb,se_parl, tvall, 2 * pnorm(abs(tvall), lower.tail = FALSE)), 5)
dimnames(ans$coefficients) <- list(names(z$coefficients),
c("Estimate", "Std. Error", "t value", "Pr(>|t|)"))
@@ -133,19 +133,19 @@
names(ans$conv_par) <- "Convergence_code_theta"
names(ans$conv_pt) <- "Sum_of_pt"
names(ans$conv_lambda) <- "Convergence_code_for_lambda"
- dimnames(ans$conv_moment) <- list(names(z$gt),"Sample_moment_with_pt")
+ dimnames(ans$conv_moment) <- list(names(z$gt), "Sample_moment_with_pt")
class(ans) <- "summary.gel"
ans
}
-residuals.gel <- function(object,...)
+residuals.gel <- function(object, ...)
{
if(is.null(object$model))
stop("The residuals method is valid only for g=formula")
object$residuals
}
-fitted.gel <- function(object,...)
+fitted.gel <- function(object, ...)
{
if(is.null(object$model))
stop("The residuals method is valid only for g=formula")
Modified: pkg/gmm/R/charStable.R
===================================================================
--- pkg/gmm/R/charStable.R 2009-12-18 01:23:31 UTC (rev 12)
+++ pkg/gmm/R/charStable.R 2009-12-21 01:10:45 UTC (rev 13)
@@ -12,7 +12,7 @@
# http://www.r-project.org/Licenses/
-charStable <- function(theta,tau,pm=0)
+charStable <- function(theta, tau, pm = 0)
{
# pm is the type parametrization as described by Nolan(2009)
# It takes the value 0 or 1
@@ -30,29 +30,29 @@
{
if(g == 0)
{
- the_car <- exp(complex(ima=d*tau))
+ the_car <- exp(complex(ima = d*tau))
}
else
{
- re_p <- -g*abs(tau)
- im_p <- d*tau
- im_p[tau!=0] <- im_p[tau!=0] + re_p[tau!=0]*2/pi*b*sign(tau[tau!=0])*log(g*abs(tau[tau!=0]))
- the_car <- exp(complex(re=re_p,ima=im_p))
+ re_p <- -g * abs(tau)
+ im_p <- d * tau
+ im_p[tau!=0] <- im_p[tau != 0] + re_p[tau != 0]*2/pi*b*sign(tau[tau != 0])*log(g*abs(tau[tau != 0]))
+ the_car <- exp(complex(re = re_p, ima = im_p))
}
}
else
{
if(g == 0)
{
- the_car <- exp(complex(ima=d*tau))
+ the_car <- exp(complex(ima = d*tau))
}
else
{
phi <- tan(pi*a/2)
re_p <- -g^a*abs(tau)^a
im_p <- d*tau*1i
- im_p[tau!=0] <- im_p[tau!=0] + re_p*( b*phi*sign(tau[tau!=0])*(abs(g*tau[tau!=0])^(1-a)-1) )
- the_car <- exp(complex(re=re_p,ima=im_p))
+ im_p[tau != 0] <- im_p[tau != 0] + re_p*( b*phi*sign(tau[tau != 0])*(abs(g*tau[tau != 0])^(1-a) - 1) )
+ the_car <- exp(complex(re = re_p, ima = im_p))
}
}
}
@@ -63,15 +63,15 @@
{
re_p <- -g*abs(tau)
im_p <- d*tau
- im_p[tau!=0] <- im_p[tau!=0]+re_p*(b*2/pi*sign(tau[tau!=0])*log(abs(tau[tau!=0])))
- the_car <- exp(complex(re=re_p,ima=im_p))
+ im_p[tau!=0] <- im_p[tau != 0]+re_p*(b*2/pi*sign(tau[tau != 0])*log(abs(tau[tau!=0])))
+ the_car <- exp(complex(re = re_p, ima = im_p))
}
else
{
phi <- tan(pi*a/2)
re_p <- -g^a*abs(tau)^a
- im_p <- re_p*(-phi*b*sign(tau))+d*tau
- the_car <- exp(complex(re=re_p,ima=im_p))
+ im_p <- re_p*(-phi*b*sign(tau)) + d*tau
+ the_car <- exp(complex(re = re_p, ima = im_p))
}
}
return(the_car)
Modified: pkg/gmm/R/gel.R
===================================================================
--- pkg/gmm/R/gel.R 2009-12-18 01:23:31 UTC (rev 12)
+++ pkg/gmm/R/gel.R 2009-12-21 01:10:45 UTC (rev 13)
@@ -11,29 +11,29 @@
# A copy of the GNU General Public License is available at
# http://www.r-project.org/Licenses/
-rho <- function(x,lamb,derive=0,type=c("EL","ET","CUE"),drop=TRUE)
+rho <- function(x, lamb, derive = 0, type = c("EL", "ET", "CUE"), drop = TRUE)
{
type <- match.arg(type)
- lamb <- matrix(lamb,ncol=1)
+ lamb <- matrix(lamb, ncol = 1)
gml <- x%*%lamb
ch <- 0
- if (derive==0)
+ if (derive == 0)
{
if (type == "EL")
{
- ch <- sum(gml>=1)
+ ch <- sum(gml >= 1)
if (drop)
{
- gml <- (gml<1)*gml
- rhomat <- log(1-gml)
+ gml <- (gml < 1)*gml
+ rhomat <- log(1 - gml)
}
else
{
- if (ch>0)
+ if (ch > 0)
rhomat <- NaN
else
- rhomat <- log(1-gml)
+ rhomat <- log(1 - gml)
}
}
if (type == "ET")
@@ -45,18 +45,18 @@
if (derive==1)
{
if (type == "EL")
- rhomat <- -1/(1-gml)
+ rhomat <- -1/(1 - gml)
if (type == "ET")
rhomat <- -exp(gml)
if (type == "CUE")
- rhomat <- -1 -gml
+ rhomat <- -1 - gml
}
if (derive==2)
{
if (type == "EL")
- rhomat <- -1/(1-gml)^2
+ rhomat <- -1/(1 - gml)^2
if (type == "ET")
rhomat <- -exp(gml)
@@ -64,14 +64,14 @@
if (type == "CUE")
rhomat <- -1
}
- rhom <-list(ch = ch, rhomat=rhomat)
+ rhom <-list(ch = ch, rhomat = rhomat)
return(rhom)
}
-getLamb <- function(g,tet,x,type=c('EL','ET','CUE'),tol_lam=1e-12,maxiterlam=1000,tol_obj=1e-7)
+getLamb <- function(g, tet, x, type = c('EL', 'ET', 'CUE'), tol_lam = 1e-12, maxiterlam = 1000, tol_obj = 1e-7)
{
type <- match.arg(type)
- gt <- g(tet,x)
+ gt <- g(tet, x)
n <- nrow(gt)
tol_cond=1e-12
gb <- colMeans(gt)
@@ -87,46 +87,46 @@
gblam0 <- NULL
j <- 1
- while ((crit > tol_lam*( 1+sqrt( crossprod(lamb0) ) ) ) & (j<=maxiterlam))
+ while ((crit > tol_lam*( 1 + sqrt( crossprod(lamb0) ) ) ) & (j <= maxiterlam))
{
- rho2 <- as.numeric(rho(gt,lamb0,derive=2,type=type)$rhomat)
- rho1 <- as.numeric(rho(gt,lamb0,derive=1,type=type)$rhomat)
+ rho2 <- as.numeric(rho(gt, lamb0, derive = 2, type = type)$rhomat)
+ rho1 <- as.numeric(rho(gt, lamb0, derive = 1, type = type)$rhomat)
gblam <- colMeans(rho1*gt)
- klam <- crossprod(rho2*gt,gt)/n
+ klam <- crossprod(rho2*gt, gt)/n
chklam <- sum(abs(klam))
if (!is.null(gblam0))
- dgblam <- crossprod(gblam)-crossprod(gblam0)
+ dgblam <- crossprod(gblam) - crossprod(gblam0)
- if (is.na(chklam) | chklam == 0 | chklam == Inf | dgblam>0 | dgblam == Inf | is.na(dgblam) | dcrit < 0)
+ if (is.na(chklam) | chklam == 0 | chklam == Inf | dgblam > 0 | dgblam == Inf | is.na(dgblam) | dcrit < 0)
{
- lamb1 <- rep(sqrt(1/n),length(lamb0))
+ lamb1 <- rep(sqrt(1/n), length(lamb0))
crit <- 0
singular=2
conv_mes <- "The algorithm produced singular system, NaN or Inf"
}
else
{
- if (rcond(klam)>tol_cond)
+ if (rcond(klam) > tol_cond)
{
- lamb1 <- lamb0-solve(klam,gblam)
- crit <- sqrt(crossprod(lamb0-lamb1))
+ lamb1 <- lamb0 - solve(klam, gblam)
+ crit <- sqrt(crossprod(lamb0 - lamb1))
lamb0 <- lamb1
}
else
{
- lamb1 <- rep(sqrt(1/n),length(lamb0))
+ lamb1 <- rep(sqrt(1/n) , length(lamb0))
crit <- 0
- singular=2
+ singular <- 2
conv_mes <- "The algorithm produced singular system"
}
}
gblam0 <- gblam
- j <- j+1
- dcrit<- crit0-crit
+ j <- j + 1
+ dcrit<- crit0 - crit
crit0 <- crit
}
- z <- list("lambda"=lamb1,singular=singular,conv_mes=conv_mes)
- if (j>maxiterlam | max(abs(gblam))>tol_obj)
+ z <- list("lambda" = lamb1, singular = singular, conv_mes = conv_mes)
+ if (j > maxiterlam | max(abs(gblam)) > tol_obj)
{
singular <- 1
conv_mes <- "No convergence after 'maxiterlam' iterations"
@@ -136,8 +136,8 @@
return(z)
}
-smoothG <- function (x, bw = bwAndrews2, prewhite = 1, ar.method = "ols",weights=weightsAndrews2,
- kernel=c("Bartlett","Parzen","Truncated","Tukey-Hanning"), approx = c("AR(1)","ARMA(1,1)"),
+smoothG <- function (x, bw = bwAndrews2, prewhite = 1, ar.method = "ols", weights = weightsAndrews2,
+ kernel = c("Bartlett", "Parzen", "Truncated", "Tukey-Hanning"), approx = c("AR(1)", "ARMA(1,1)"),
tol = 1e-7)
{
kernel <- match.arg(kernel)
@@ -157,31 +157,31 @@
rt <- length(w)
if (rt >= 2)
{
- w <- c(w[rt:2],w)
+ w <- c(w[rt:2], w)
w <- w / sum(w)
- rt <- rt-1
- sgt <- function(t) crossprod(x[(t-rt):(t+rt),],w)
- x[(rt+1):(n-rt),] <- t(sapply((rt+1):(n-rt),sgt))
- sx <- list("smoothx"=x,"kern_weights"=w)
+ rt <- rt - 1
+ sgt <- function(t) crossprod(x[(t-rt):(t+rt),], w)
+ x[(rt+1):(n-rt),] <- t(sapply((rt + 1):(n - rt), sgt))
+ sx <- list("smoothx" = x, "kern_weights" = w)
return(sx)
}
else
- sx <- list("smoothx"=x,"kern_weights"=1)
+ sx <- list("smoothx" = x,"kern_weights" = 1)
return(sx)
}
-gel <- function(g,x,tet0,gradv=NULL,smooth=FALSE,type=c("EL","ET","CUE","ETEL"), vcov=c("HAC","iid"), kernel = c("Bartlett", "Parzen",
- "Truncated", "Tukey-Hanning"), bw=bwAndrews2, approx = c("AR(1)",
- "ARMA(1,1)"), prewhite = 1, ar.method = "ols", tol_weights = 1e-7, tol_lam=1e-9, tol_obj = 1e-9,
- tol_mom = 1e-9,maxiterlam=1000, constraint=FALSE,optfct=c("optim","optimize","nlminb"),optlam=c("iter","numeric"),
- model=TRUE, X=FALSE, Y=FALSE,...)
+gel <- function(g, x, tet0, gradv = NULL, smooth = FALSE, type = c("EL", "ET", "CUE", "ETEL"), vcov = c("HAC", "iid"),
+ kernel = c("Bartlett", "Parzen", "Truncated", "Tukey-Hanning"), bw = bwAndrews2, approx = c("AR(1)",
+ "ARMA(1,1)"), prewhite = 1, ar.method = "ols", tol_weights = 1e-7, tol_lam = 1e-9, tol_obj = 1e-9,
+ tol_mom = 1e-9, maxiterlam = 1000, constraint = FALSE, optfct = c("optim", "optimize", "nlminb"),
+ optlam = c("iter", "numeric"), model = TRUE, X = FALSE, Y = FALSE, ...)
{
- vcov=match.arg(vcov)
+ vcov <- match.arg(vcov)
type <- match.arg(type)
optfct <- match.arg(optfct)
optlam <- match.arg(optlam)
- weights=weightsAndrews2
+ weights <- weightsAndrews2
if (type == "ETEL")
{
typel <- "ET"
@@ -194,57 +194,58 @@
}
approx <- match.arg(approx)
kernel <- match.arg(kernel)
- if(optfct=="optim")
+ if(optfct == "optim")
k <- length(tet0)
else
k <- 1
typeg=0
- if (is(g,"formula"))
+ if (is(g, "formula"))
{
- typeg=1
- dat <- getDat(g,x)
+ typeg = 1
+ dat <- getDat(g, x)
x <- dat$x
- g <- function(tet,x,ny=dat$ny,nh=dat$nh,k=dat$k)
+ g <- function(tet, x, ny = dat$ny, nh = dat$nh, k = dat$k)
{
- tet <- matrix(tet,ncol=k)
- e <- x[,1:ny] - x[,(ny+1):(ny+k)]%*%t(tet)
- gt <- e*x[,ny+k+1]
+ tet <- matrix(tet, ncol = k)
+ e <- x[,1:ny] - x[, (ny+1):(ny+k)]%*%t(tet)
+ gt <- e*x[, ny+k+1]
if (nh > 1)
{
for (i in 2:nh)
{
- gt <- cbind(gt,e*x[,(ny+k+i)])
+ gt <- cbind(gt, e*x[,(ny+k+i)])
}
}
return(gt)
}
- gradv <- function(tet,x,ny=dat$ny,nh=dat$nh,k=dat$k)
+ gradv <- function(tet, x, ny = dat$ny, nh = dat$nh, k = dat$k)
{
- tet <- matrix(tet,ncol=k)
- dgb <- -(t(x[,(ny+k+1):(ny+k+nh)])%*%x[,(ny+1):(ny+k)])%x%diag(rep(1,ny))/nrow(x)
+ tet <- matrix(tet, ncol = k)
+ dgb <- -(t(x[,(ny+k+1):(ny+k+nh)])%*%x[,(ny+1):(ny+k)])%x%diag(rep(1, ny))/nrow(x)
return(dgb)
}
}
if (typeg)
n <- nrow(x)
else
- n = nrow(g(tet0,x))
+ n = nrow(g(tet0, x))
if (smooth)
{
g1 <- g
- rgmm <- gmm(g,x,tet0,wmatrix="ident")
+ rgmm <- gmm(g, x, tet0, wmatrix = "ident")
if (is.function(weights))
- w <- weights(g(rgmm$coefficients,x),kernel=kernel,bw=bw,prewhite = prewhite,ar.method=ar.method,approx=approx,tol=tol_weights)
+ w <- weights(g(rgmm$coefficients, x), kernel = kernel, bw = bw, prewhite = prewhite,
+ ar.method = ar.method, approx = approx, tol = tol_weights)
else
w <- weights
- sg <- function(thet,x)
+ sg <- function(thet, x)
{
- gf <- g1(thet,x)
- gt <- smoothG(gf, weights=w)$smoothx
+ gf <- g1(thet, x)
+ gt <- smoothG(gf, weights = w)$smoothx
return(gt)
}
g <- sg
@@ -252,41 +253,41 @@
lll <- 1
thetf <- function(tet)
{
- if (optlam=="iter")
+ if (optlam == "iter")
{
- lamblist <- getLamb(g,tet,x,type=typel,tol_lam=tol_lam,maxiterlam=maxiterlam,tol_obj=tol_obj)
+ lamblist <- getLamb(g, tet, x, type = typel, tol_lam = tol_lam, maxiterlam = maxiterlam, tol_obj = tol_obj)
lamb <- lamblist$lambda
- gt <- g(tet,x)
- pt <- -rho(gt,lamb,type=typet,derive=1)$rhomat/nrow(gt)
+ gt <- g(tet, x)
+ pt <- -rho(gt, lamb, type = typet, derive = 1)$rhomat/nrow(gt)
checkmom <- sum(as.numeric(pt)*gt)
- if (lamblist$singular==0)
- p <- sum(rho(gt,lamb,type=typet)$rhomat) + abs(checkmom)/tol_mom
- if (lamblist$singular==1)
- p <- sum(rho(gt,lamb,type=typet)$rhomat) + abs(checkmom)/tol_mom + lamblist$obj/tol_mom
- if (lamblist$singular==2)
+ if (lamblist$singular == 0)
+ p <- sum(rho(gt, lamb, type = typet)$rhomat) + abs(checkmom)/tol_mom
+ if (lamblist$singular == 1)
+ p <- sum(rho(gt, lamb, type = typet)$rhomat) + abs(checkmom)/tol_mom + lamblist$obj/tol_mom
+ if (lamblist$singular == 2)
p <- 1e50*lll
- lll <- lll+1
+ lll <- lll + 1
}
else
{
- gt <- g(tet,x)
+ gt <- g(tet, x)
rhofct <- function(lamb)
{
- rhof <- -sum(rho(gt,lamb,type=typel)$rhomat)
+ rhof <- -sum(rho(gt, lamb, type = typel)$rhomat)
return(rhof)
}
- if (ncol(gt)>1)
- rlamb <- optim(rep(0,ncol(gt)),rhofct,control=list(maxit=1000))
+ if (ncol(gt) > 1)
+ rlamb <- optim(rep(0, ncol(gt)), rhofct, control = list(maxit = 1000))
else
{
- rlamb <- optimize(rhofct,c(-1,1))
+ rlamb <- optimize(rhofct, c(-1,1))
rlamb$par <- rlamb$minimum
rlamb$value <- rlamb$objective
}
lamb <- rlamb$par
- pt <- -rho(gt,lamb,type=typet,derive=1)$rhomat/nrow(gt)
+ pt <- -rho(gt, lamb, type = typet, derive = 1)$rhomat/nrow(gt)
checkmom <- sum(as.numeric(pt)*gt)
- p <- -rlamb$value + (checkmom)^2/tol_mom + (sum(as.numeric(pt))-1)^2/tol_mom
+ p <- -rlamb$value + (checkmom)^2/tol_mom + (sum(as.numeric(pt)) - 1)^2/tol_mom
}
return(p)
}
@@ -294,49 +295,49 @@
if (!constraint)
{
if (optfct == "optim")
- res <- optim(tet0,thetf,...)
+ res <- optim(tet0, thetf, ...)
if (optfct == "nlminb")
- res <- nlminb(tet0,thetf,...)
+ res <- nlminb(tet0, thetf, ...)
if (optfct == "optimize")
{
- res <- optimize(thetf,tet0,...)
+ res <- optimize(thetf, tet0, ...)
res$par <- res$minimum
res$convergence <- "There is no convergence code for optimize"
}
}
if(constraint)
- res <- constrOptim(tet0,thetf,grad=NULL,...)
+ res <- constrOptim(tet0, thetf, grad = NULL, ...)
if (optlam=="iter")
{
- rlamb <- getLamb(g,res$par,x,type=typel,tol_lam=tol_lam,maxiterlam=maxiterlam,tol_obj=tol_obj)
- z <- list(coefficients=res$par,lambda=rlamb$lam,conv_lambda=rlamb$conv_mes,conv_par=res$convergence)
+ rlamb <- getLamb(g,res$par, x, type = typel, tol_lam = tol_lam, maxiterlam = maxiterlam, tol_obj = tol_obj)
+ z <- list(coefficients = res$par, lambda = rlamb$lam, conv_lambda = rlamb$conv_mes, conv_par = res$convergence)
z$foc_lambda <- rlamb$obj
}
if (optlam=="numeric")
{
- gt<-g(res$par,x)
+ gt<-g(res$par, x)
rhofct <- function(lamb)
{
- rhof <- -sum(rho(gt,lamb,type=typel)$rhomat)
+ rhof <- -sum(rho(gt, lamb, type = typel)$rhomat)
return(rhof)
}
- rlamb <- optim(rep(0,ncol(gt)),rhofct,control=list(maxit=1000))
- z <- list(coefficients=res$par,conv_par=res$convergence,lambda=rlamb$par)
- z$conv_lambda=paste("Lambda by optim. Conv. code = ",rlamb$convergence,sep="")
- rho1 <- as.numeric(rho(gt,z$lambda,derive=1,type=typel)$rhomat)
+ rlamb <- optim(rep(0, ncol(gt)), rhofct, control = list(maxit = 1000))
+ z <- list(coefficients = res$par, conv_par = res$convergence, lambda = rlamb$par)
+ z$conv_lambda = paste("Lambda by optim. Conv. code = ", rlamb$convergence, sep = "")
+ rho1 <- as.numeric(rho(gt, z$lambda, derive = 1, type = typel)$rhomat)
z$foc_lambda <- crossprod(colMeans(rho1*gt))
}
z$type <- type
- z$gt <- g(z$coefficients,x)
- rhom <- rho(z$gt,z$lambda,type=typet)
- z$pt <- -rho(z$gt,z$lambda,type=typet,derive=1)$rhomat/n
+ z$gt <- g(z$coefficients, x)
+ rhom <- rho(z$gt, z$lambda, type = typet)
+ z$pt <- -rho(z$gt, z$lambda, type = typet, derive = 1)$rhomat/n
z$conv_moment <- colSums(as.numeric(z$pt)*z$gt)
z$conv_pt <- sum(as.numeric(z$pt))
- z$objective <- sum(as.numeric(rhom$rhomat)-rho(1,0,type=typet)$rhomat)/n
+ z$objective <- sum(as.numeric(rhom$rhomat) - rho(1, 0, type = typet)$rhomat)/n
if (type == "EL")
{
@@ -347,34 +348,35 @@
if (!is.function(gradv))
G <- .Gf(z$coefficients, x, g)
else
- G <- gradv(z$coefficients,x)
+ G <- gradv(z$coefficients, x)
if (vcov == "iid")
khat <- crossprod(z$gt)
else
- khat <- HAC(g(z$coefficients,x), kernel=kernel, bw=bw,prewhite=prewhite,ar.method=ar.method,approx=approx,tol=tol_weights)
+ khat <- HAC(g(z$coefficients, x), kernel = kernel, bw = bw, prewhite = prewhite,
+ ar.method=ar.method, approx = approx, tol = tol_weights)
- kg <- solve(khat,G)
- z$vcov_par <- solve(crossprod(G,kg))/n
- z$vcov_lambda <- ((solve(khat)-kg%*%z$vcov_par%*%t(kg)))/n
+ kg <- solve(khat, G)
+ z$vcov_par <- solve(crossprod(G, kg))/n
+ z$vcov_lambda <- ((solve(khat) - kg%*%z$vcov_par%*%t(kg)))/n
- if (smooth) z$weights<-w
+ if (smooth) z$weights <- w
if (typeg ==0)
{
- names(z$coefficients) <- paste("Theta[",1:k,"]",sep="")
- colnames(z$gt) <- paste("gt[",1:ncol(z$gt),"]",sep="")
- names(z$lambda) <- paste("Lambda[",1:ncol(z$gt),"]",sep="")
+ names(z$coefficients) <- paste("Theta[",1:k,"]", sep = "")
+ colnames(z$gt) <- paste("gt[",1:ncol(z$gt),"]", sep = "")
+ names(z$lambda) <- paste("Lambda[",1:ncol(z$gt),"]", sep = "")
}
if (typeg == 1)
{
- namex <- colnames(dat$x[,(dat$ny+1):(dat$ny+dat$k)])
- nameh <- colnames(dat$x[,(dat$ny+dat$k+1):(dat$ny+dat$k+dat$nh)])
+ namex <- colnames(dat$x[, (dat$ny+1):(dat$ny+dat$k)])
+ nameh <- colnames(dat$x[, (dat$ny+dat$k+1):(dat$ny+dat$k+dat$nh)])
if (dat$ny > 1)
{
- namey <- colnames(dat$x[,1:dat$ny])
- names(z$coefficients) <- paste(rep(namey,dat$k),"_",rep(namex,rep(dat$ny,dat$k)),sep="")
- colnames(z$gt) <- paste(rep(namey,dat$nh),"_",rep(nameh,rep(dat$ny,dat$nh)),sep="")
- names(z$lambda) <- paste("Lam(",rep(namey,dat$nh),"_",rep(nameh,rep(dat$ny,dat$nh)),")",sep="")
+ namey <- colnames(dat$x[, 1:dat$ny])
+ names(z$coefficients) <- paste(rep(namey, dat$k), "_", rep(namex, rep(dat$ny, dat$k)), sep = "")
+ colnames(z$gt) <- paste(rep(namey, dat$nh), "_", rep(nameh, rep(dat$ny,dat$nh)), sep = "")
+ names(z$lambda) <- paste("Lam(",rep(namey,dat$nh), "_", rep(nameh, rep(dat$ny,dat$nh)), ")", sep = "")
}
if (dat$ny == 1)
{
@@ -383,25 +385,25 @@
names(z$lambda) <- nameh
}
}
- dimnames(z$vcov_par) <- list(names(z$coefficients),names(z$coefficients))
- dimnames(z$vcov_lambda) <- list(names(z$lambda),names(z$lambda))
- if (typeg==1)
+ dimnames(z$vcov_par) <- list(names(z$coefficients), names(z$coefficients))
+ dimnames(z$vcov_lambda) <- list(names(z$lambda), names(z$lambda))
+ if (typeg == 1)
{
b <- z$coefficients
y <- as.matrix(model.response(dat$mf, "numeric"))
ny <- dat$ny
- b <- t(matrix(b,nrow=ny))
+ b <- t(matrix(b, nrow = ny))
x <- as.matrix(model.matrix(dat$mt, dat$mf, NULL))
yhat <- x%*%b
- z$fitted.values<-yhat
- z$residuals<-y-yhat
- z$terms<- dat$mt
- if(model) z$model<-dat$mf
- if(X) z$x<-x
- if(Y) z$y<-y
+ z$fitted.values <- yhat
+ z$residuals <- y - yhat
+ z$terms <- dat$mt
+ if(model) z$model <- dat$mf
+ if(X) z$x <- x
+ if(Y) z$ y<- y
}
else
- if(X) z$x<-x
+ if(X) z$x <- x
z$call <- match.call()
class(z) <- "gel"
Modified: pkg/gmm/R/specTest.R
===================================================================
--- pkg/gmm/R/specTest.R 2009-12-18 01:23:31 UTC (rev 12)
+++ pkg/gmm/R/specTest.R 2009-12-21 01:10:45 UTC (rev 13)
@@ -7,9 +7,9 @@
{
j <- x$objective*x$n
J_test <- noquote(paste("J-Test: degrees of freedom is ",x$df,sep=""))
- j <- noquote(cbind(j,ifelse(x$df>0,pchisq(j,x$df,lower.tail = FALSE),"*******")))
- dimnames(j) <- list("Test E(g)=0: ",c("J-test","P-value"))
- ans<-list(ntest=J_test,test=j)
+ j <- noquote(cbind(j, ifelse(x$df>0,pchisq(j,x$df, lower.tail = FALSE),"*******")))
+ dimnames(j) <- list("Test E(g)=0: ", c("J-test", "P-value"))
+ ans<-list(ntest=J_test, test = j)
class(ans) <- "specTest"
ans
}
@@ -29,15 +29,15 @@
khat <- crossprod(x$gt)/n
gbar <- colMeans(x$gt)
LR_test <- 2*x$objective*n
- LM_test <- n*crossprod(x$lambda,crossprod(khat,x$lambda))
- J_test <- n*crossprod(gbar,solve(khat,gbar))
- test <- c(LR_test,LM_test,J_test)
- df <- (ncol(x$gt)-length(x$par))
- ntest <- noquote(paste("Over-identifying restrictions tests: degrees of freedom is ",df,sep=""))
- vptest <- pchisq(test,df,lower.tail=FALSE)
+ LM_test <- n*crossprod(x$lambda, crossprod(khat, x$lambda))
+ J_test <- n*crossprod(gbar, solve(khat, gbar))
+ test <- c(LR_test, LM_test, J_test)
+ df <- (ncol(x$gt) - length(x$par))
+ ntest <- noquote(paste("Over-identifying restrictions tests: degrees of freedom is ", df, sep = ""))
+ vptest <- pchisq(test,df,lower.tail = FALSE)
test <- cbind(test,vptest)
- dimnames(test) <- list(c("LR test","LM test","J test"),c("statistics","p-value"))
- ans <- list(test=test,ntest=ntest)
+ dimnames(test) <- list(c("LR test", "LM test", "J test"), c("statistics", "p-value"))
+ ans <- list(test = test, ntest = ntest)
class(ans) <- "specTest"
ans
}
More information about the Gmm-commits
mailing list