1# This program is free software; you can redistribute it and/or modify 2# it under the terms of the GNU General Public License as published by 3# the Free Software Foundation; either version 2 of the License, or 4# (at your option) any later version. 5# 6# This program is distributed in the hope that it will be useful, 7# but WITHOUT ANY WARRANTY; without even the implied warranty of 8# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 9# GNU General Public License for more details. 10# 11# A copy of the GNU General Public License is available at 12# http://www.r-project.org/Licenses/ 13 14gmm <- function(g,x,t0=NULL,gradv=NULL, type=c("twoStep","cue","iterative"), wmatrix = c("optimal","ident"), vcov=c("HAC","MDS","iid","TrueFixed"), 15 kernel=c("Quadratic Spectral","Truncated", "Bartlett", "Parzen", "Tukey-Hanning"),crit=10e-7,bw = bwAndrews, 16 prewhite = 1, ar.method = "ols", approx="AR(1)",tol = 1e-7, itermax=100,optfct=c("optim","optimize","nlminb", "constrOptim"), 17 model=TRUE, X=FALSE, Y=FALSE, TypeGmm = "baseGmm", centeredVcov = TRUE, weightsMatrix = NULL, traceIter = FALSE, data, eqConst = NULL, 18 eqConstFullVcov = FALSE, mustar = NULL, onlyCoefficients=FALSE, ...) 19 { 20 type <- match.arg(type) 21 kernel <- match.arg(kernel) 22 vcov <- match.arg(vcov) 23 wmatrix <- match.arg(wmatrix) 24 optfct <- match.arg(optfct) 25 26 if (!is.null(eqConst)) 27 TypeGmm <- "constGmm" 28 29 if(vcov=="TrueFixed" & is.null(weightsMatrix)) 30 stop("TrueFixed vcov only for fixed weighting matrix") 31 if(!is.null(weightsMatrix)) 32 wmatrix <- "optimal" 33 34 if(missing(data)) 35 data<-NULL 36 all_args<-list(data = data, g = g, x = x, t0 = t0, gradv = gradv, type = type, wmatrix = wmatrix, vcov = vcov, kernel = kernel, 37 crit = crit, bw = bw, prewhite = prewhite, ar.method = ar.method, approx = approx, 38 weightsMatrix = weightsMatrix, centeredVcov = centeredVcov, tol = tol, itermax = itermax, 39 optfct = optfct, model = model, X = X, Y = Y, call = match.call(), traceIter = traceIter, 40 eqConst = eqConst, eqConstFullVcov = eqConstFullVcov, mustar = mustar) 41 class(all_args)<-TypeGmm 42 Model_info<-getModel(all_args, ...) 43 z <- momentEstim(Model_info, ...) 44 if (onlyCoefficients) 45 return(z[c("coefficients","objective")]) 46 z <- FinRes(z, Model_info) 47 z 48 } 49 50evalGmm <- function(g, x, t0, tetw=NULL, gradv=NULL, wmatrix = c("optimal","ident"), 51 vcov=c("HAC","iid","TrueFixed"), 52 kernel=c("Quadratic Spectral","Truncated", "Bartlett", "Parzen", 53 "Tukey-Hanning"),crit=10e-7,bw = bwAndrews, 54 prewhite = FALSE, ar.method = "ols", approx="AR(1)",tol = 1e-7, 55 model=TRUE, X=FALSE, Y=FALSE, centeredVcov = TRUE, 56 weightsMatrix = NULL, data, mustar = NULL) 57 { 58 TypeGmm = "baseGmm" 59 type <- "eval" 60 kernel <- match.arg(kernel) 61 vcov <- match.arg(vcov) 62 wmatrix <- match.arg(wmatrix) 63 if (is.null(tetw) & is.null(weightsMatrix)) 64 stop("If the weighting matrix is not provided, you need to provide the vector of parameters tetw") 65 if(vcov=="TrueFixed" & is.null(weightsMatrix)) 66 stop("TrueFixed vcov only for fixed weighting matrix") 67 if(!is.null(weightsMatrix)) 68 wmatrix <- "optimal" 69 if(missing(data)) 70 data<-NULL 71 all_args<-list(data = data, g = g, x = x, t0 = t0, tetw = tetw, gradv = gradv, 72 type = type, wmatrix = wmatrix, vcov = vcov, kernel = kernel, 73 crit = crit, bw = bw, prewhite = prewhite, ar.method = ar.method, 74 approx = approx, weightsMatrix = weightsMatrix, 75 centeredVcov = centeredVcov, tol = tol, itermax = 100, 76 model = model, X = X, Y = Y, call = match.call(), 77 traceIter = NULL, optfct="optim", 78 eqConst = NULL, eqConstFullVcov = FALSE, mustar = mustar) 79 class(all_args)<-TypeGmm 80 Model_info<-getModel(all_args) 81 class(Model_info) <- "baseGmm.eval" 82 z <- momentEstim(Model_info) 83 z <- FinRes(z, Model_info) 84 z 85 } 86 87tsls <- function(g,x,data) 88 { 89 if(class(g)[1] != "formula") 90 stop("2SLS is for linear models expressed as formula only") 91 ans <- gmm(g,x,data=data,vcov="iid", TypeGmm="tsls") 92 ans$met <- "Two Stage Least Squares" 93 ans$call <- match.call() 94 class(ans) <- c("tsls","gmm") 95 ans$vcov <- vcov(ans, type="Classical") 96 return(ans) 97 } 98 99.myKernHAC <- function(gmat, obj) 100 { 101 gmat <- as.matrix(gmat) 102 if(obj$centeredVcov) 103 gmat <- scale(gmat, scale=FALSE) 104 class(gmat) <- "gmmFct" 105 AllArg <- obj$WSpec$sandwich 106 AllArg$x <- gmat 107 if (is.function(AllArg$bw)) 108 { 109 if (identical(AllArg$bw, bwWilhelm)) 110 { 111 G <- .DmomentFct(obj$tet, obj$dat) 112 obj <- list(gt=gmat, G=G) 113 class(obj) <- "gmm" 114 } else { 115 obj <- gmat 116 } 117 AllArg$bw <- AllArg$bw(obj, order.by = AllArg$order.by, 118 kernel = AllArg$kernel, 119 prewhite = AllArg$prewhite, 120 ar.method = AllArg$ar.method, 121 approx=AllArg$approx) 122 } 123 weights <- weightsAndrews(x=gmat, bw=AllArg$bw, kernel=AllArg$kernel, 124 prewhite=AllArg$prewhite, tol=AllArg$tol) 125 w <- vcovHAC(x=gmat, order.by=AllArg$order.by, weights=weights, 126 prewhite=AllArg$prewhite, sandwich=FALSE, 127 ar.method=AllArg$ar.method, adjust=FALSE) 128 attr(w,"Spec") <- list(weights = weights, bw = AllArg$bw, 129 kernel = AllArg$kernel) 130 w 131 } 132 133getDat <- function (formula, h, data, error=TRUE) 134{ 135 cl <- match.call() 136 mf <- match.call(expand.dots = FALSE) 137 m <- match(c("formula", "data"), names(mf), 0L) 138 mf <- mf[c(1L, m)] 139 mf$drop.unused.levels <- TRUE 140 mf$na.action <- "na.pass" 141 mf[[1L]] <- as.name("model.frame") 142 mf <- eval(mf, parent.frame()) 143 mt <- attr(mf, "terms") 144 y <- as.matrix(model.response(mf, "numeric")) 145 namey <- as.character(formula)[2] 146 if (ncol(y)>1) 147 namey <- paste(namey, ".", 1:ncol(y), sep="") 148 xt <- as.matrix(model.matrix(mt, mf, NULL)) 149 n <- NROW(y) 150 if (inherits(h,'formula')) 151 { 152 tmp <- as.character(formula) 153 termsh <- terms(h) 154 h <- paste(tmp[2], "~", as.character(h)[2], sep="") 155 h <- as.formula(h) 156 mfh <- match.call(expand.dots = FALSE) 157 mh <- match(c("h", "data"), names(mfh), 0L) 158 mfh <- mfh[c(1L, mh)] 159 mfh$formula <- h 160 mfh$h <- NULL 161 mfh$drop.unused.levels <- TRUE 162 mfh$na.action <- "na.pass" 163 mfh[[1L]] <- as.name("model.frame") 164 mfh <- eval(mfh, parent.frame()) 165 mth <- attr(mfh, "terms") 166 h <- as.matrix(model.matrix(mth, mfh, NULL)) 167 } 168 else 169 { 170 h <- as.matrix(h) 171 chkInt <- sapply(1:ncol(h), function(i) all(h[,i]/mean(h[,i]) == 1)) 172 if (sum(chkInt) > 1) 173 stop("Too many intercept in the matrix h") 174 if (any(chkInt)) 175 { 176 h <- h[,!chkInt, drop=FALSE] 177 h <- cbind(1,h) 178 intercept_h <- TRUE 179 } else { 180 if (attr(mt,"intercept")==1) 181 { 182 h <- cbind(1, h) 183 intercept_h <- TRUE 184 } else { 185 intercept_h <- FALSE 186 } 187 } 188 if(is.null(colnames(h))) 189 { 190 if (intercept_h) 191 colnames(h) <- c("(Intercept)",paste("h",1:(ncol(h)-1),sep="")) 192 else 193 colnames(h) <- paste("h",1:ncol(h),sep="") 194 } else { 195 if (intercept_h) 196 colnames(h)[1] <- "(Intercept)" 197 coln_h <- colnames(h) 198 coln_h <- gsub(" ", "", coln_h) 199 chk <- which(coln_h == "") 200 if (length(chk) >0) 201 coln_h[chk] <- paste("h", 1:length(chk), sep="") 202 if (any(duplicated(coln_h))) 203 stop("colnames of the matrix h must be unique") 204 colnames(h) <- coln_h 205 } 206 if (!intercept_h) 207 { 208 hFormula <- paste(colnames(h), collapse="+") 209 hFormula <- as.formula(paste("~", hFormula, "-1", sep="")) 210 } else { 211 hFormula <- paste(colnames(h)[-1], collapse="+") 212 hFormula <- as.formula(paste("~", hFormula, sep="")) 213 } 214 termsh <- terms(hFormula) 215 } 216 ny <- ncol(y) 217 k <- ncol(xt) 218 nh <- ncol(h) 219 if (nh<k) 220 { 221 if (error) 222 stop("The number of moment conditions must be at least equal to the number of coefficients to estimate") 223 } 224 if (is.null(colnames(y))) 225 colnames(y) <- namey 226 rownames(xt) <- rownames(y) 227 rownames(h) <- rownames(y) 228 x <- cbind(y,xt,h) 229 if(any(is.na(x))) 230 { 231 warning("There are missing values. Associated observations have been removed") 232 x <- na.omit(x) 233 if (nrow(x)<=k) 234 { 235 if (error) 236 stop("The number of observations must be greater than the number of coefficients") 237 } 238 } 239 colnames(x)<-c(colnames(y),colnames(xt),colnames(h)) 240 return(list(x=x,nh=nh,ny=ny,k=k,mf=mf,mt=mt,cl=cl,termsh=termsh,termsx=mt)) 241} 242 243.tetlin <- function(dat, w, type=NULL) 244 { 245 x <- dat$x 246 g <- .momentFct 247 gradv <- .DmomentFct 248 ny <- dat$ny 249 nh <- dat$nh 250 k <- dat$k 251 n <- nrow(x) 252 ym <- as.matrix(x[,1:ny]) 253 xm <- as.matrix(x[,(ny+1):(ny+k)]) 254 hm <- as.matrix(x[,(ny+k+1):(ny+k+nh)]) 255 if (!is.null(attr(dat, "eqConst"))) 256 { 257 resTet <- attr(dat,"eqConst")$eqConst 258 y2 <- xm[, resTet[,1],drop=FALSE]%*%resTet[,2] 259 ym <- ym-c(y2) 260 xm <- xm[,-resTet[,1],drop=FALSE] 261 k <- ncol(xm) 262 } 263 includeExo <- which(colnames(xm)%in%colnames(hm)) 264 inv <- attr(w, "inv") 265 mustar <- attr(dat, "mustar") 266 if (!is.null(type)) 267 { 268 if(type=="2sls") 269 { 270 if (length(includeExo) > 0) 271 { 272 endo <- xm[, -includeExo, drop = FALSE] 273 endoName <- colnames(endo) 274 if (ncol(endo) != 0) 275 { 276 if (attr(dat$termsh, "intercept") == 1) 277 restsls <- lm(endo~hm[,-1]) 278 else 279 restsls <- lm(endo~hm-1) 280 fsls <- xm 281 fsls[, -includeExo] <- restsls$fitted 282 } else { 283 fsls <- xm 284 restsls <- NULL 285 } 286 } else { 287 if (attr(dat$termsh, "intercept") == 1) 288 restsls <- lm(xm~hm[,-1]) 289 else 290 restsls <- lm(xm~hm-1) 291 fsls <- restsls$fitted 292 endoName <- colnames(xm) 293 } 294 par <- lm.fit(as.matrix(fsls), ym)$coefficients 295 if (ny == 1) 296 { 297 e2sls <- ym-xm%*%par 298 v2sls <- lm.fit(as.matrix(hm), e2sls)$fitted 299 value <- sum(v2sls^2)/sum(e2sls^2) 300 } else { 301 par <- c(t(par)) 302 g2sls <- g(par, dat) 303 w <- crossprod(g2sls)/n 304 gb <- matrix(colMeans(g2sls), ncol = 1) 305 value <- crossprod(gb, solve(w, gb)) 306 } 307 } 308 } else { 309 if (ny>1) 310 { 311 if (inv) 312 { 313 whx <- solve(w, (crossprod(hm, xm) %x% diag(ny))) 314 wvecyh <- solve(w, matrix(crossprod(ym, hm), ncol = 1)) 315 } else { 316 whx <- w%*% (crossprod(hm, xm) %x% diag(ny)) 317 wvecyh <- w%*%matrix(crossprod(ym, hm), ncol = 1) 318 } 319 dg <- gradv(NULL, dat) 320 xx <- crossprod(dg, whx) 321 par <- solve(xx, crossprod(dg, wvecyh)) 322 } else { 323 if (nh>k) 324 { 325 if(inv) 326 xzwz <- crossprod(xm,hm)%*%solve(w,t(hm)) 327 else 328 xzwz <- crossprod(xm,hm)%*%w%*%t(hm) 329 par <- solve(xzwz%*%xm,xzwz%*%ym) 330 } else { 331 par <- solve(crossprod(hm,xm),crossprod(hm,ym)) 332 } 333 } 334 gb <- matrix(colSums(g(par, dat))/n, ncol = 1) 335 if(inv) 336 value <- crossprod(gb, solve(w, gb)) 337 else 338 value <- crossprod(gb, w%*%gb) 339 } 340 res <- list(par = par, value = value) 341 if (!is.null(mustar)) 342 { 343 if (!is.null(type)) 344 { 345 w <- crossprod(hm)/NROW(hm) 346 attr(w, "inv") <- TRUE 347 } 348 res <- .mustarLin(par, xm, hm, w, dat, mustar) 349 } 350 if (!is.null(type)) 351 { 352 if (type == "2sls") 353 res$firstStageReg <- restsls 354 if (!is.null(restsls)) 355 { 356 chk <- .chkPerfectFit(restsls) 357 res$fsRes <- suppressWarnings(summary(restsls))[!chk] 358 attr(res$fsRes, "Endo") <- endoName[!chk] 359 } 360 } 361 return(res) 362 } 363 364.mustarLin <- function(par, xm, hm, w, dat, mustar) 365 { 366 if (ncol(xm) == ncol(hm)) 367 { 368 par <- par-solve(crossprod(hm,xm),mustar) 369 } else { 370 hmxm <- crossprod(hm,xm) 371 if (attr(w, "inv")) 372 T1 <- solve(w, hmxm) 373 else 374 T1 <- w%*%hmxm 375 par <- par-solve(crossprod(hmxm, T1), crossprod(T1, mustar)) 376 } 377 gb <- colSums(.momentFct(par, dat))/NROW(xm) 378 if(attr(w, "inv")) 379 value <- crossprod(gb, solve(w, gb)) 380 else 381 value <- crossprod(gb, w%*%gb) 382 list(value=value, par=par) 383 } 384 385.obj1 <- function(thet, x, w) 386 { 387 gt <- .momentFct(thet, x) 388 gbar <- as.vector(colMeans(gt)) 389 INV <- attr(w, "inv") 390 if (INV) 391 obj <- crossprod(gbar, solve(w, gbar)) 392 else 393 obj <- crossprod(gbar,w)%*%gbar 394 return(obj) 395 } 396 397.Gf <- function(thet, x, pt = NULL) 398 { 399 myenv <- new.env() 400 assign('x', x, envir = myenv) 401 assign('thet', thet, envir = myenv) 402 barg <- function(x, thet) 403 { 404 gt <- .momentFct(thet, x) 405 if (is.null(pt)) 406 gbar <- as.vector(colMeans(gt)) 407 else 408 gbar <- as.vector(colSums(c(pt)*gt)) 409 410 return(gbar) 411 } 412 G <- attr(numericDeriv(quote(barg(x, thet)), "thet", myenv), "gradient") 413 return(G) 414 } 415 416.objCue <- function(thet, x, type=c("HAC", "MDS", "iid", "ident", "fct", "fixed")) 417 { 418 type <- match.arg(type) 419 gt <- .momentFct(thet, x) 420 gbar <- as.vector(colMeans(gt)) 421 w <- .weightFct(thet, x, type) 422 if (attr(w, "inv")) 423 obj <- crossprod(gbar,solve(w,gbar)) 424 else 425 obj <- crossprod(gbar,w%*%gbar) 426 return(obj) 427 } 428 429 430.momentFct <- function(tet, dat) 431 { 432 if (!is.null(attr(dat, "eqConst"))) 433 { 434 resTet <- attr(dat,"eqConst")$eqConst 435 tet2 <- vector(length=length(tet)+nrow(resTet)) 436 tet2[resTet[,1]] <- resTet[,2] 437 tet2[-resTet[,1]] <- tet 438 } else { 439 tet2 <- tet 440 } 441 if (attr(dat, "ModelType") == "linear") 442 { 443 x <- dat$x 444 ny <- dat$ny 445 nh <- dat$nh 446 k <- dat$k 447 tet2 <- matrix(tet2, ncol = k) 448 e <- x[,1:ny] - x[,(ny+1):(ny+k)] %*% t(tet2) 449 gt <- e * x[, ny+k+1] 450 if(nh > 1) 451 for (i in 2:nh) 452 gt <- cbind(gt, e*x[, (ny+k+i)]) 453 } else { 454 gt <- attr(dat,"momentfct")(tet2, dat) 455 } 456 if (!is.null(attr(dat, "smooth"))) 457 { 458 bw <- attr(dat, "smooth")$bw 459 w <- attr(dat, "smooth")$w 460 gt <- smoothG(gt, bw = bw, weights = w)$smoothx 461 } 462 gt <- as.matrix(gt) 463 if (!is.null(attr(dat, "mustar"))) 464 { 465 if (length(attr(dat, "mustar")) != ncol(gt)) 466 stop("dim of mustar must match the number of moment conditions") 467 gt <- sweep(gt, 2, attr(dat, "mustar"), "-") 468 } 469 return(gt) 470 } 471 472.DmomentFct <- function(tet, dat, pt=NULL) 473 { 474 if (!is.null(attr(dat, "eqConst"))) 475 { 476 resTet <- attr(dat,"eqConst")$eqConst 477 tet2 <- vector(length=length(tet)+nrow(resTet)) 478 tet2[resTet[,1]] <- resTet[,2] 479 tet2[-resTet[,1]] <- tet 480 } else { 481 tet2 <- tet 482 } 483 if ((attr(dat,"ModelType") == "linear") & (is.null(attr(dat, "smooth")))) 484 { 485 x <- dat$x 486 ny <- dat$ny 487 nh <- dat$nh 488 k <- dat$k 489 dgb <- -(t(x[,(ny+k+1):(ny+k+nh)]) %*% x[,(ny+1):(ny+k)]) %x% diag(rep(1,ny))/nrow(x) 490 if (!is.null(attr(dat, "eqConst"))) 491 dgb <- dgb[,-attr(dat,"eqConst")$eqConst[,1], drop=FALSE] 492 } else { 493 if (is.null(attr(dat,"gradv"))) 494 { 495 dgb <- .Gf(tet, dat, pt) 496 } else { 497 dgb <- attr(dat,"gradv")(tet2, dat) 498 if (!is.null(attr(dat, "eqConst"))) 499 dgb <- dgb[,-attr(dat,"eqConst")$eqConst[,1], drop=FALSE] 500 } 501 } 502 return(as.matrix(dgb)) 503 } 504 505.weightFct <- function(tet, dat, type=c("HAC", "MDS", "iid", "ident", "fct", "fixed")) 506 { 507 type <- match.arg(type) 508 if (type == "fixed") 509 { 510 w <- attr(dat, "weight")$w 511 attr(w, "inv") <- FALSE 512 } else if (type == "ident") { 513 w <- diag(attr(dat, "q")) 514 attr(w, "inv") <- FALSE 515 } else { 516 gt <- .momentFct(tet,dat) 517 if (!is.null(attr(dat, "namesgt"))) 518 colnames(gt) <- attr(dat, "namesgt") 519 if(attr(dat, "weight")$centeredVcov) 520 gt <- scale(gt, scale=FALSE) 521 n <- NROW(gt) 522 } 523 if (type == "HAC") 524 { 525 obj <- attr(dat, "weight") 526 obj$centeredVcov <- FALSE 527 obj$tet <- tet 528 obj$dat <- dat 529 w <- .myKernHAC(gt, obj) 530 attr(w, "inv") <- TRUE 531 } 532 if (type == "MDS") 533 { 534 w <- crossprod(gt)/n 535 attr(w, "inv") <- TRUE 536 } 537 if (type == "iid") 538 { 539 if (attr(dat, "ModelType") == "linear") 540 { 541 if (dat$ny == 1) 542 { 543 e <- .residuals(tet, dat)$residuals 544 sig <- mean(scale(e,scale=FALSE)^2) 545 z <- dat$x[,(1+dat$ny+dat$k):ncol(dat$x)] 546 w <- sig*crossprod(z)/length(e) 547 } else { 548 w <- crossprod(gt)/n 549 } 550 } else { 551 w <- crossprod(gt)/n 552 } 553 attr(w, "inv") <- TRUE 554 } 555 if (type == "fct") 556 { 557 w <- attr(dat, "weight")$fct(gt, attr(dat, "weight")$fctArg) 558 attr(w, "inv") <- TRUE 559 } 560 return(w) 561 } 562 563.residuals <- function(tet, dat) 564 { 565 if (!is.null(attr(dat, "eqConst"))) 566 { 567 resTet <- attr(dat,"eqConst")$eqConst 568 tet2 <- vector(length=length(tet)+nrow(resTet)) 569 tet2[resTet[,1]] <- resTet[,2] 570 tet2[-resTet[,1]] <- tet 571 } else { 572 tet2 <- tet 573 } 574 tet2 <- t(matrix(tet2, nrow = dat$ny)) 575 y <- as.matrix(dat$x[,1:dat$ny]) 576 x <- as.matrix(dat$x[,(dat$ny+1):(dat$ny+dat$k)]) 577 yhat <- x%*%tet2 578 e <- y-yhat 579 return(list(residuals=e, yhat=yhat)) 580 } 581 582