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 14.rho <- function(x, lamb, derive = 0, type = c("EL", "ET", "CUE", "HD", "ETEL", "ETHD", 15 "RCUE"), k = 1) 16 { 17 type <- match.arg(type) 18 if (type == "RCUE") 19 type <- "CUE" 20 gml <- x%*%c(lamb)*k 21 if (derive == 0) 22 { 23 if (type == "EL") 24 { 25 if (any(gml>=1)) 26 stop("Computation of Lambda fails because NAs produced by log(1-gt*l)") 27 rhomat <- log(1 - gml) 28 } 29 if (type == "ET") 30 rhomat <- -exp(gml) 31 if (type == "CUE") 32 rhomat <- -gml -0.5*gml^2 33 if (type == "HD") 34 rhomat <- -2/(1-gml/2) 35 if (type == "ETEL") 36 { 37 w <- -exp(gml) 38 w <- w/sum(w) 39 rhomat <- -log(w*NROW(gml)) 40 } 41 if (type == "ETHD") 42 { 43 w <- -exp(gml) 44 w <- w/sum(w) 45 rhomat <- (sqrt(w)-1/sqrt(NROW(gml)))^2 46 } 47 } 48 if (derive==1) 49 { 50 if (type == "EL") 51 rhomat <- -1/(1 - gml) 52 if (type == "ET") 53 rhomat <- -exp(gml) 54 if (type == "CUE") 55 rhomat <- -1 - gml 56 if (type == "HD") 57 rhomat <- -1/((1-gml/2)^2) 58 if (any(type == c("ETEL","ETHD"))) 59 rhomat <- NULL 60 } 61 if (derive==2) 62 { 63 if (type == "EL") 64 rhomat <- -1/(1 - gml)^2 65 if (type == "ET") 66 rhomat <- -exp(gml) 67 if (type == "CUE") 68 rhomat <- -rep(1,nrow(x)) 69 if (type == "HD") 70 rhomat <- -1/((1-gml/2)^3) 71 if (any(type == c("ETEL","ETHD"))) 72 rhomat <- NULL 73 } 74 return(c(rhomat)) 75 } 76 77.getCgelLam <- function(gt, l0, type = c('EL', 'ET', 'CUE', "HD"), 78 method = c("nlminb", "optim", "constrOptim"), 79 control=list(), k = 1, alpha = 0.01) 80 { 81 type <- match.arg(type) 82 method <- match.arg(method) 83 fct <- function(l, X) 84 { 85 r1 <- colMeans(.rho(gt,l,derive=1,type=type,k=k)*X) 86 crossprod(r1) + alpha*crossprod(l) 87 } 88 Dfct <- function(l, X) 89 { 90 r2 <- .rho(X,l,derive=2,type=type,k=k) 91 r1 <- .rho(X,l,derive=1,type=type,k=k) 92 H <- t(X*r2)%*%X/nrow(X) 93 2*H%*%colMeans(r1*X) + 2*alpha*l 94 } 95 96 if (method == "nlminb") 97 res <- nlminb(l0, fct, Dfct, X = gt, control = control) 98 if (method == "optim") 99 res <- optim(l0, fct, Dfct, X = gt, method="BFGS", control = control) 100 if (method == "constrOptim") 101 { 102 ci <- -rep(1,nrow(gt)) 103 res <- constrOptim(rep(0,ncol(gt)),fct,Dfct,-gt,ci,control=control,X=gt) 104 } 105 if (method == "optim") 106 { 107 conv <- list(convergence = res$convergence, 108 counts = res$counts, message = res$message) 109 } else { 110 conv <- list(convergence = res$convergence, counts = res$evaluations, 111 message = res$message) 112 } 113 114 return(list(lambda = res$par, convergence = conv, 115 obj = mean(.rho(gt,res$par, derive=0,type=type,k=k)) - 116 .rho(1,0, derive=0,type=type,k=k))) 117 } 118 119.Wu <- function(gt, tol_lam = 1e-8, maxiterlam = 50, K=1) 120 { 121 u <- as.matrix(gt) 122 n=length(nrow(u)) 123 M=rep(0,ncol(u)) 124 dif=1 125 tol=tol_lam 126 k=0 127 while(dif>tol & k<=maxiterlam) 128 { 129 D1=t(u)%*%((1/(1+u%*%M))*rep(1,n)) 130 DD=-t(u)%*%(c((1/(1+u%*%M)^2))*u) 131 D2=solve(DD,D1,tol=1e-40) 132 dif=max(abs(D2)) 133 rule=1 134 while(rule>0) 135 { 136 rule=0 137 if(min(1+t(M-D2)%*%t(u))<=0) rule=rule+1 138 if(rule>0) D2=D2/2 139 } 140 M=M-D2 141 k=k+1 142 } 143 if(k>=maxiterlam) 144 { 145 M=rep(0,ncol(u)) 146 conv = list(convergence=1) 147 } else { 148 conv = list(convergence=0) 149 } 150 return(list(lambda = c(-M), convergence = conv, obj = 151 mean(.rho(gt,-M,derive=0,type="EL",k=K)))) 152 } 153 154.Wu2 <- function(gt, tol_lam = 1e-8, maxiter = 50, K = 1) 155 { 156 gt <- as.matrix(gt) 157 res <- .Fortran(F_wu, as.double(gt), as.double(tol_lam), 158 as.integer(maxiter), as.integer(nrow(gt)), 159 as.integer(ncol(gt)), as.double(K), 160 conv=integer(1), obj=double(1), 161 lambda=double(ncol(gt))) 162 list(lambda=res$lambda, convergence=list(convergence=res$conv), 163 obj = res$obj) 164 } 165 166.CUE_lam <- function(gt, K=1) 167 { 168 q <- qr(gt) 169 n <- nrow(gt) 170 l0 <- -qr.coef(q, rep(1,n)) 171 conv <- list(convergence=0) 172 list(lambda = l0, convergence = conv, obj = 173 mean(.rho(gt,l0,derive=0,type="CUE",k=K))) 174 } 175 176.CUE_lamPos <- function(gt, K=1, control=list()) 177 { 178 getpt <- function(gt,lam) 179 { 180 gl <- c(gt%*%lam) 181 pt <- 1 + gl 182 pt/sum(pt) 183 } 184 maxit <- ifelse("maxit" %in% names(control), 185 control$maxit, 50) 186 res <-.CUE_lam(gt, K) 187 n <- nrow(gt) 188 i <- 1 189 pt <- getpt(gt, res$lambda) 190 w <- pt<0 191 while (TRUE) 192 { 193 gt2 <- gt[!w,] 194 n1 <- nrow(gt2) 195 if (n1 == n) 196 break 197 res <- try(.CUE_lam(gt2), silent=TRUE) 198 if (i > maxit) 199 return(list(lambda=rep(0,ncol(gt)), obj=0, pt=rep(1/n,n), 200 convergence=list(convergence=1))) 201 if (any(class(res) == "try-error")) 202 return(list(lambda=rep(0,ncol(gt)), obj=0, pt=rep(1/n,n), 203 convergence=list(convergence=2))) 204 pt[!w] <- getpt(gt2, res$lambda) 205 pt[w] <- 0 206 if (all(pt>=0)) 207 break 208 i <- i+1 209 w[!w] <- pt[!w]<0 210 } 211 n0 <- n-n1 212 res$obj <- res$obj*n1/n + n0/(2*n) 213 res 214 } 215 216.CUE_lamPos2 <- function(gt, K=1, control=list()) 217 { 218 gt <- as.matrix(gt) 219 n <- nrow(gt) 220 q <- ncol(gt) 221 maxit <- ifelse("maxit" %in% names(control), 222 control$maxit, 50) 223 res <- try(.Fortran(F_lamcuep, as.double(gt), 224 as.integer(n), as.integer(q), as.double(K), 225 as.integer(maxit),conv=integer(1), 226 lam=double(q),pt=double(n), 227 obj=double(1) 228 ), silent=TRUE) 229 if (any(class(res) == "try-error")) 230 return(list(lambda=rep(0,q), obj=0, pt=rep(1/n,n), 231 convergence=list(convergence=3))) 232 list(lambda=res$lam, obj=res$obj, pt=res$pt, 233 convergence=list(convergence=res$conv)) 234 } 235 236getLamb <- function(gt, l0, type = c('EL', 'ET', 'CUE', "ETEL", "HD", "ETHD", "RCUE"), 237 tol_lam = 1e-7, maxiterlam = 100, tol_obj = 1e-7, 238 k = 1, method = c("nlminb", "optim", "iter", "Wu"), 239 control=list()) 240 { 241 method <- match.arg(method) 242 type <- match.arg(type) 243 gt <- as.matrix(gt) 244 if (method == "Wu" & type != "EL") 245 stop("Wu (2005) method to compute Lambda is for EL only") 246 if (method == "Wu") 247 return(.Wu2(gt, tol_lam, maxiterlam, k)) 248 if (type == "CUE") 249 return(.CUE_lam(gt, k)) 250 if (type == "RCUE") 251 return(.CUE_lamPos2(gt, k, control)) 252 if (method == "iter") 253 { 254 if ((type == "ETEL") | (type == "ETHD")) 255 type = "ET" 256 for (i in 1:maxiterlam) 257 { 258 r1 <- .rho(gt,l0,derive=1,type=type,k=k) 259 r2 <- .rho(gt,l0,derive=2,type=type,k=k) 260 F <- -colMeans(r1*gt) 261 J <- crossprod(r2*gt,gt) 262 if (sum(abs(F))<tol_obj) 263 { 264 conv <- list(convergence="Tolerance for the FOC reached") 265 break 266 } 267 P <- solve(J,F) 268 if (sum(abs(P))<tol_lam) 269 { 270 conv <- list(convergence="Tolerance on lambda reached") 271 break 272 } 273 l0 <- l0 + P 274 conv <- list(convergence="maxiterlam reached") 275 } 276 } else { 277 fct <- function(l,X,type,k) 278 { 279 r0 <- .rho(X,l,derive=0,type=type,k=k) 280 -mean(r0) 281 } 282 Dfct <- function(l,X,type,k) 283 { 284 r1 <- .rho(X,l,derive=1,type=type,k=k) 285 -colMeans(r1*X) 286 } 287 DDfct <- function(l,X,type,k) 288 { 289 r2 <- .rho(X,l,derive=2,type=type,k=k) 290 -t(X*r2)%*%X/nrow(X) 291 } 292 if ((type == "ETEL")|(type=="ETHD")) 293 type = "ET" 294 if (method=="optim") 295 { 296 if (type != "EL") 297 { 298 res <- optim(rep(0,ncol(gt)),fct,gr=Dfct,X=gt,type=type, 299 k=k,method="BFGS",control=control) 300 } else { 301 ci <- -rep(1,nrow(gt)) 302 res <- constrOptim(rep(0,ncol(gt)),fct,Dfct,-gt,ci, 303 control=control,X=gt,type=type,k=k) 304 } 305 } else { 306 res <- nlminb(rep(0,ncol(gt)), fct, gradient = Dfct, 307 hessian = DDfct, X = gt, type=type, k=k, 308 control = control) 309 } 310 l0 <- res$par 311 if (method == "optim" | method == "constrOptim") 312 conv <- list(convergence = res$convergence, counts = res$counts, 313 message = res$message) 314 if(method == "nlminb") 315 conv <- list(convergence = res$convergence, counts = 316 res$evaluations, message = res$message) 317 } 318 return(list(lambda = l0, convergence = conv, obj = 319 mean(.rho(gt,l0,derive=0,type=type,k=k))- 320 .rho(1, 0, type = type, derive = 0, k = k))) 321 } 322 323smoothG <- function (x, bw = bwAndrews, prewhite = 1, ar.method = "ols", 324 weights = weightsAndrews, 325 kernel = c("Bartlett", "Parzen", "Truncated", "Tukey-Hanning"), 326 approx = c("AR(1)", "ARMA(1,1)"), 327 tol = 1e-7) 328 { 329 kernel <- match.arg(kernel) 330 approx <- match.arg(approx) 331 332 n <- nrow(x) 333 if (is.function(weights)) 334 { 335 class(x) <- "gmmFct" 336 w <- weights(x, bw = bw, kernel = kernel, 337 prewhite = prewhite, ar.method = ar.method, tol = tol, 338 verbose = FALSE, approx = approx) 339 } else { 340 w <- weights 341 } 342 if (is.numeric(w)) 343 { 344 rt <- length(w) 345 if (rt >= 2) 346 { 347 w <- c(w[rt:2], w) 348 w <- w / sum(w) 349 w <- kernel(w[rt:length(w)]) 350 } else { 351 w <- kernel(1) 352 } 353 } else { 354 if (class(w)[1] != "tskernel") 355 stop("Provided weights must be a numeric vector or an object of class 'tskernel'") 356 } 357 if (length(w$coef)>1) 358 x <- kernapply(x, w) 359 sx <- list("smoothx" = x, "kern_weights" = w, bw = bw) 360 return(sx) 361 } 362 363gel <- function(g, x, tet0 = NULL, gradv = NULL, smooth = FALSE, 364 type = c("EL", "ET", "CUE", "ETEL", "HD", "ETHD","RCUE"), 365 kernel = c("Truncated", "Bartlett"), bw = bwAndrews, 366 approx = c("AR(1)", "ARMA(1,1)"), prewhite = 1, ar.method = "ols", 367 tol_weights = 1e-7, tol_lam = 1e-9, tol_obj = 1e-9, 368 tol_mom = 1e-9, maxiterlam = 100, constraint = FALSE, 369 optfct = c("optim", "optimize", "nlminb"), 370 optlam = c("nlminb", "optim", "iter", "Wu"), data, 371 Lambdacontrol = list(), 372 model = TRUE, X = FALSE, Y = FALSE, TypeGel = "baseGel", alpha = NULL, 373 eqConst = NULL, eqConstFullVcov = FALSE, onlyCoefficients=FALSE, ...) 374 { 375 type <- match.arg(type) 376 optfct <- match.arg(optfct) 377 optlam <- match.arg(optlam) 378 weights <- weightsAndrews 379 approx <- match.arg(approx) 380 kernel <- match.arg(kernel) 381 if (!is.null(eqConst)) 382 TypeGel <- "constGel" 383 384 if(missing(data)) 385 data<-NULL 386 all_args <- list(g = g, x = x, tet0 = tet0, gradv = gradv, smooth = smooth, 387 type = type, kernel = kernel, bw = bw, approx = approx, 388 prewhite = prewhite, ar.method = ar.method, 389 tol_weights = tol_weights, tol_lam = tol_lam, 390 tol_obj = tol_obj, tol_mom = tol_mom, maxiterlam = maxiterlam, 391 constraint = constraint, optfct = optfct, weights = weights, 392 optlam = optlam, model = model, X = X, Y = Y, 393 TypeGel = TypeGel, call = match.call(), 394 Lambdacontrol = Lambdacontrol, alpha = alpha, data = data, 395 eqConst = eqConst, eqConstFullVcov = eqConstFullVcov, 396 onlyCoefficients=onlyCoefficients) 397 class(all_args)<-TypeGel 398 Model_info<-getModel(all_args) 399 z <- momentEstim(Model_info, ...) 400 if (!onlyCoefficients) 401 class(z) <- "gel" 402 return(z) 403 } 404 405evalGel <- function(g, x, tet0, gradv = NULL, smooth = FALSE, 406 type = c("EL", "ET", "CUE", "ETEL", "HD", "ETHD","RCUE"), 407 kernel = c("Truncated", "Bartlett"), bw = bwAndrews, 408 approx = c("AR(1)", "ARMA(1,1)"), prewhite = 1, 409 ar.method = "ols", tol_weights = 1e-7, tol_lam = 1e-9, 410 tol_obj = 1e-9, tol_mom = 1e-9, maxiterlam = 100, 411 optlam = c("nlminb", "optim", "iter", "Wu"), data, 412 Lambdacontrol = list(), 413 model = TRUE, X = FALSE, Y = FALSE, alpha = NULL, ...) 414 { 415 type <- match.arg(type) 416 optlam <- match.arg(optlam) 417 weights <- weightsAndrews 418 approx <- match.arg(approx) 419 kernel <- match.arg(kernel) 420 TypeGel <- "baseGel" 421 422 if(missing(data)) 423 data<-NULL 424 all_args <- list(g = g, x = x, tet0 = tet0, gradv = gradv, smooth = smooth, 425 type = type, kernel = kernel, bw = bw, approx = approx, 426 prewhite = prewhite, ar.method = ar.method, 427 tol_weights = tol_weights, tol_lam = tol_lam, 428 tol_obj = tol_obj, tol_mom = tol_mom, maxiterlam = maxiterlam, 429 weights = weights, optlam = optlam, model = model, X = X, 430 Y = Y, call = match.call(), Lambdacontrol = Lambdacontrol, 431 alpha = alpha, data = data, optfct="optim") 432 class(all_args)<-TypeGel 433 Model_info<-getModel(all_args) 434 class(Model_info) <- "baseGel.eval" 435 z <- momentEstim(Model_info, ...) 436 class(z) <- "gel" 437 return(z) 438 } 439 440.thetf <- function(tet, P, output=c("obj","all"), l0Env) 441 { 442 output <- match.arg(output) 443 gt <- P$g(tet, P$x) 444 l0 <- get("l0",envir=l0Env) 445 if (((P$type=="ETEL")|(P$type=="ETHD"))&(!is.null(P$CGEL))) 446 { 447 P$CGEL <- NULL 448 warning("CGEL not implemented for ETEL or for ETHD") 449 } 450 if (is.null(P$CGEL)) 451 { 452 if (P$optlam != "optim" & P$type == "EL") 453 { 454 lamb <- try(getLamb(gt, l0, type = P$type, tol_lam = P$tol_lam, 455 maxiterlam = P$maxiterlam, 456 tol_obj = P$tol_obj, k = P$k1/P$k2, 457 control = P$Lambdacontrol, 458 method = P$optlam), silent = TRUE) 459 if(any(class(lamb) == "try-error")) 460 lamb <- getLamb(gt, l0, type = P$type, tol_lam = P$tol_lam, 461 maxiterlam = P$maxiterlam, 462 tol_obj = P$tol_obj, k = P$k1/P$k2, 463 control = P$Lambdacontrol, method = "optim") 464 } else { 465 lamb <- getLamb(gt, l0, type = P$type, tol_lam = P$tol_lam, 466 maxiterlam = P$maxiterlam, tol_obj = P$tol_obj, 467 k = P$k1/P$k2, control = P$Lambdacontrol, 468 method = P$optlam) 469 } 470 } else { 471 lamb <- try(.getCgelLam(gt, l0, type = P$type, method = "nlminb", 472 control=P$Lambdacontrol, k = P$k1/P$k2, 473 alpha = P$CGEL),silent=TRUE) 474 if (any(class(lamb) == "try-error")) 475 lamb <- try(.getCgelLam(gt, l0, type = P$type, 476 method = "constrOptim", 477 control=P$Lambdacontrol, 478 k = P$k1/P$k2, alpha = P$CGEL),silent=TRUE) 479 } 480 if (P$type == "ETEL") 481 obj <- mean(.rho(gt, lamb$lambda, type = P$type, derive = 0, 482 k = P$k1/P$k2) - .rho(1, 0, type = P$type, 483 derive = 0, k = P$k1/P$k2)) 484 else if (P$type == "ETHD") 485 obj <- sum(.rho(gt, lamb$lambda, type = P$type, derive = 0, 486 k = P$k1/P$k2) - .rho(1, 0, type = P$type, derive = 0, 487 k = P$k1/P$k2)) 488 else 489 obj <- lamb$obj 490 assign("l0",lamb$lambda,envir=l0Env) 491 if(output == "obj") 492 return(obj) 493 else 494 return(list(obj = obj, lambda = lamb, gt = gt)) 495 } 496 497.getImpProb <- function(gt, lambda, type, k1, k2) 498 { 499 if ((type == "ETEL")|(type=="ETHD")) 500 type <- "ET" 501 n <- NROW(gt) 502 pt <- -.rho(gt, lambda, type = type, derive = 1, k = k1/k2)/n 503 # Making sure pt>0 504 if (type=="CUE") 505 { 506 eps <- -length(pt)*min(min(pt),0) 507 pt <- (pt+eps/length(pt))/(1+eps) 508 } 509 if (type=="RCUE") 510 pt[pt<0] <- 0 511 ################### 512 conv_moment <- colSums(pt*gt) 513 conv_pt <- sum(as.numeric(pt)) 514 pt <- pt/sum(pt) 515 attr(pt, "conv_moment") <- conv_moment 516 attr(pt, "conv_pt") <- conv_pt 517 pt 518 } 519 520.vcovGel <- function(gt, G, k1, k2, bw, pt=NULL,tol=1e-16) 521 { 522 q <- NCOL(gt) 523 n <- NROW(gt) 524 if (is.null(pt)) 525 pt <- 1/n 526 G <- G/k1 527 gt <- gt*sqrt(pt*bw/k2) 528 qrGt <- qr(gt) 529 piv <- sort.int(qrGt$pivot, index.return=TRUE)$ix 530 R <- qr.R(qrGt)[,piv] 531 X <- forwardsolve(t(R), G) 532 Y <- forwardsolve(t(R), diag(q)) 533 res <- lm.fit(X,Y) 534 u <- res$residuals 535 Sigma <- chol2inv(res$qr$qr)/n 536 diag(Sigma)[diag(Sigma)<0] <- tol 537 if (q==ncol(G)) 538 { 539 SigmaLam <- matrix(0, q, q) 540 } else { 541 SigmaLam <- backsolve(R, u)/n*bw^2 542 diag(SigmaLam)[diag(SigmaLam)<0] <- tol 543 } 544 khat <- crossprod(R) 545 list(vcov_par=Sigma, vcov_lambda=SigmaLam,khat=khat) 546 } 547