1# These functions are 2# Copyright (C) 1998-2021 T.W. Yee, University of Auckland. 3# All rights reserved. 4 5 6 7 8 9 10 11 12 13rrar.Ci <- function(i, coeffs, aa, Ranks., MM) { 14 index <- cumsum(c(aa, MM*Ranks.)) 15 ans <- matrix(coeffs[(index[i]+1):index[i+1]], 16 Ranks.[i], MM, byrow = TRUE) 17 t(ans) 18} 19 20 21rrar.Ak1 <- function(MM, coeffs, Ranks., aa) { 22 ptr <- 0 23 Ak1 <- diag(MM) 24 for (jay in 1:MM) { 25 for (i in 1:MM) { 26 if (i > jay && (MM+1)-(Ranks.[jay]-1) <= i) { 27 ptr <- ptr + 1 28 Ak1[i,jay] <- coeffs[ptr] 29 } 30 } 31 } 32 if (aa > 0 && ptr != aa) stop("something wrong") 33 Ak1 34} 35 36 37rrar.Di <- function(i, Ranks.) { 38 if (Ranks.[1] == Ranks.[i]) diag(Ranks.[i]) else 39 rbind(diag(Ranks.[i]), 40 matrix(0, Ranks.[1] - Ranks.[i], Ranks.[i])) 41} 42 43 44rrar.Mi <- function(i, MM, Ranks., ki) { 45 if (Ranks.[ki[i]] == MM) 46 return(NULL) 47 hi <- Ranks.[ki[i]] - Ranks.[ki[i+1]] 48 Ji <- matrix(0, hi, Ranks.[1]) 49 for (j in 1:hi) { 50 Ji[j,j+Ranks.[ki[i+1]]] <- 1 51 } 52 Mi <- matrix(0, MM-Ranks.[ki[i]], MM) # dim(Oi) == dim(Ji) 53 for (j in 1:(MM-Ranks.[ki[i]])) { 54 Mi[j,j+Ranks.[ki[i ]]] <- 1 55 } 56 kronecker(Mi, Ji) 57} 58 59 60rrar.Mmat <- function(MM, uu, Ranks., ki) { 61 Mmat <- NULL 62 for (ii in uu:1) { 63 Mmat <- rbind(Mmat, rrar.Mi(ii, MM, Ranks., ki)) 64 } 65 Mmat 66} 67 68 69block.diag <- function(A, B) { 70 if (is.null(A) && is.null(B)) 71 return(NULL) 72 if (!is.null(A) && is.null(B)) 73 return(A) 74 if (is.null(A) && !is.null(B)) 75 return(B) 76 77 A <- as.matrix(A) 78 B <- as.matrix(B) 79 temp <- cbind(A, matrix(0, nrow(A), ncol(B))) 80 rbind(temp, cbind(matrix(0, nrow(B), ncol(A)), B)) 81} 82 83 84rrar.Ht <- function(plag, MM, Ranks., coeffs, aa, uu, ki) { 85 Htop <- Hbot <- NULL 86 Mmat <- rrar.Mmat(MM, uu, Ranks., ki) # NULL if full rank 87 Ak1 <- rrar.Ak1(MM, coeffs, Ranks., aa) 88 89 if (!is.null(Mmat)) 90 for (i in 1:plag) { 91 Di <- rrar.Di(i, Ranks.) 92 Ci <- rrar.Ci(i, coeffs, aa, Ranks., MM) 93 temp <- Di %*% t(Ci) 94 Htop <- cbind(Htop, Mmat %*% kronecker(diag(MM), temp)) 95 } 96 97 for (i in 1:plag) { 98 Di <- rrar.Di(i, Ranks.) 99 temp <- kronecker(t(Di) %*% t(Ak1), diag(MM)) 100 Hbot <- block.diag(Hbot, temp) 101 } 102 rbind(Htop, Hbot) 103} 104 105 106rrar.Ut <- function(y, tt, plag, MM) { 107 Ut <- NULL 108 if (plag>1) 109 for (i in 1:plag) { 110 Ut <- rbind(Ut, kronecker(diag(MM), cbind(y[tt-i,]))) 111 } 112 Ut 113} 114 115 116rrar.UU <- function(y, plag, MM, n) { 117 UU <- NULL 118 for (i in (plag+1):n) { 119 UU <- rbind(UU, t(rrar.Ut(y, i, plag, MM))) 120 } 121 UU 122} 123 124 125rrar.Wmat <- function(y, Ranks., MM, ki, plag, aa, uu, n, coeffs) { 126 temp1 <- rrar.UU(y, plag, MM, n) 127 temp2 <- t(rrar.Ht(plag, MM, Ranks., coeffs, aa, uu, ki)) 128 list(UU = temp1, Ht = temp2) 129} 130 131 132 133 134 135rrar.control <- 136 function(stepsize = 0.5, save.weights = TRUE, 137 summary.HDEtest = FALSE, # Overwrites the summary() default. 138 ...) { 139 140 if (stepsize <= 0 || stepsize > 1) { 141 warning("bad value of stepsize; using 0.5 instead") 142 stepsize <- 0.5 143 } 144 list(stepsize = stepsize, 145 summary.HDEtest = summary.HDEtest, 146 save.weights = as.logical(save.weights)[1]) 147} 148 149 150 151 rrar <- function(Ranks = 1, coefstart = NULL) { 152 lag.p <- length(Ranks) 153 154 new("vglmff", 155 blurb = c("Nested reduced-rank vector autoregressive model AR(", 156 lag.p, ")\n\n", 157 "Link: ", 158 namesof("mu_t", "identitylink"), 159 ", t = ", paste(paste(1:lag.p, coll = ",", sep = ""))), 160 initialize = eval(substitute(expression({ 161 Ranks. <- .Ranks 162 plag <- length(Ranks.) 163 nn <- nrow(x) # original n 164 indices <- 1:plag 165 166 copy.X.vlm <- TRUE # X.vlm.save matrix changes at each iteration 167 168 dsrank <- -sort(-Ranks.) # ==rev(sort(Ranks.)) 169 if (any(dsrank != Ranks.)) 170 stop("Ranks must be a non-increasing sequence") 171 if (!is.matrix(y) || ncol(y) == 1) { 172 stop("response must be a matrix with more than one column") 173 } else { 174 MM <- ncol(y) 175 ki <- udsrank <- unique(dsrank) 176 uu <- length(udsrank) 177 for (i in 1:uu) 178 ki[i] <- max((1:plag)[dsrank == udsrank[i]]) 179 ki <- c(ki, plag+1) # For computing a 180 Ranks. <- c(Ranks., 0) # For computing a 181 aa <- sum( (MM-Ranks.[ki[1:uu]]) * 182 (Ranks.[ki[1:uu]]-Ranks.[ki[-1]]) ) 183 } 184 if (!intercept.only) 185 warning("ignoring explanatory variables") 186 187 if (any(MM < Ranks.)) 188 stop("'max(Ranks)' can only be ", MM, " or less") 189 y.save <- y # Save the original 190 if (any(w != 1)) 191 stop("all weights should be 1") 192 193 new.coeffs <- .coefstart # Needed for iter = 1 of $weight 194 new.coeffs <- if (length(new.coeffs)) 195 rep_len(new.coeffs, aa+sum(Ranks.)*MM) else 196 runif(aa+sum(Ranks.)*MM) 197 temp8 <- rrar.Wmat(y.save, Ranks., MM, ki, plag, 198 aa, uu, nn, new.coeffs) 199 X.vlm.save <- temp8$UU %*% temp8$Ht 200 201 if (!length(etastart)) { 202 etastart <- X.vlm.save %*% new.coeffs 203 etastart <- matrix(etastart, ncol = ncol(y), byrow = TRUE) 204 } 205 206 extra$Ranks. <- Ranks.; extra$aa <- aa 207 extra$plag <- plag; extra$nn <- nn 208 extra$MM <- MM; extra$coeffs <- new.coeffs; 209 extra$y.save <- y.save 210 211 keep.assign <- attr(x, "assign") 212 x <- x[-indices, , drop = FALSE] 213 if (is.R()) 214 attr(x, "assign") <- keep.assign 215 y <- y[-indices, , drop = FALSE] 216 w <- w[-indices] 217 n.save <- n <- nn - plag 218 }), list( .Ranks = Ranks, .coefstart = coefstart ))), 219 220 linkinv = function(eta, extra = NULL) { 221 aa <- extra$aa 222 coeffs <- extra$coeffs 223 MM <- extra$MM 224 nn <- extra$nn 225 plag <- extra$plag 226 Ranks. <- extra$Ranks. 227 y.save <- extra$y.save 228 229 tt <- (1+plag):nn 230 mu <- matrix(0, nn-plag, MM) 231 Ak1 <- rrar.Ak1(MM, coeffs, Ranks., aa) 232 for (i in 1:plag) { 233 Di <- rrar.Di(i, Ranks.) 234 Ci <- rrar.Ci(i, coeffs, aa, Ranks., MM) 235 mu <- mu + y.save[tt-i, , drop = FALSE] %*% 236 t(Ak1 %*% Di %*% t(Ci)) 237 } 238 mu 239 }, 240 last = expression({ 241 misc$plag <- plag 242 misc$Ranks <- Ranks. 243 misc$Ak1 <- Ak1 244 misc$omegahat <- omegahat 245 misc$Cmatrices <- Cmatrices 246 misc$Dmatrices <- Dmatrices 247 misc$Hmatrix <- temp8$Ht 248 misc$Phimatrices <- vector("list", plag) 249 for (ii in 1:plag) { 250 misc$Phimatrices[[ii]] <- Ak1 %*% Dmatrices[[ii]] %*% 251 t(Cmatrices[[ii]]) 252 } 253 misc$Z <- y.save %*% t(solve(Ak1)) 254 }), 255 vfamily = "rrar", 256 validparams = function(eta, y, extra = NULL) { 257 okay1 <- TRUE 258 okay1 259 }, 260 deriv = expression({ 261 temp8 <- rrar.Wmat(y.save, Ranks., MM, ki, plag, 262 aa, uu, nn, new.coeffs) 263 X.vlm.save <- temp8$UU %*% temp8$Ht 264 265 extra$coeffs <- new.coeffs 266 267 resmat <- y 268 tt <- (1+plag):nn 269 Ak1 <- rrar.Ak1(MM, new.coeffs, Ranks., aa) 270 Cmatrices <- Dmatrices <- vector("list", plag) 271 for (ii in 1:plag) { 272 Dmatrices[[ii]] <- Di <- rrar.Di(ii, Ranks.) 273 Cmatrices[[ii]] <- Ci <- rrar.Ci(ii, new.coeffs, aa, Ranks., MM) 274 resmat <- resmat - y.save[tt - ii, , drop = FALSE] %*% 275 t(Ak1 %*% Di %*% t(Ci)) 276 } 277 omegahat <- (t(resmat) %*% resmat) / n # MM x MM 278 omegainv <- solve(omegahat) 279 280 omegainv <- solve(omegahat) 281 ind1 <- iam(NA, NA, MM, both = TRUE) 282 283 wz <- matrix(omegainv[cbind(ind1$row, ind1$col)], 284 nn-plag, length(ind1$row), byrow = TRUE) 285 mux22(t(wz), y-mu, M = extra$MM, as.matrix = TRUE) 286 }), 287 weight = expression({ 288 wz 289 })) 290} 291 292 293 294 295 296 297 298 299vglm.garma.control <- function(save.weights = TRUE, ...) { 300 list(save.weights = as.logical(save.weights)[1]) 301} 302 303 304 garma <- function(link = "identitylink", 305 p.ar.lag = 1, 306 q.ma.lag = 0, 307 coefstart = NULL, 308 step = 1.0) { 309 310 if (!is.Numeric(p.ar.lag, integer.valued = TRUE, length.arg = 1)) 311 stop("bad input for argument 'p.ar.lag'") 312 if (!is.Numeric(q.ma.lag, integer.valued = TRUE, length.arg = 1)) 313 stop("bad input for argument 'q.ma.lag'") 314 if (q.ma.lag != 0) 315 stop("sorry, only q.ma.lag = 0 is currently implemented") 316 317 318 link <- as.list(substitute(link)) 319 earg <- link2list(link) 320 link <- attr(earg, "function.name") 321 322 323 new("vglmff", 324 blurb = c("GARMA(", p.ar.lag, ",", q.ma.lag, ")\n\n", 325 "Link: ", 326 namesof("mu_t", link, earg = earg), 327 ", t = ", paste(paste(1:p.ar.lag, coll = ",", sep = ""))), 328 initialize = eval(substitute(expression({ 329 plag <- .p.ar.lag 330 predictors.names <- namesof("mu", .link , earg = .earg , tag = FALSE) 331 332 indices <- 1:plag 333 tt.index <- (1 + plag):nrow(x) 334 p.lm <- ncol(x) 335 336 copy.X.vlm <- TRUE # x matrix changes at each iteration 337 338 if ( .link == "logitlink" || .link == "probitlink" || 339 .link == "clogloglink" || .link == "cauchitlink") { 340 delete.zero.colns <- TRUE 341 eval(process.categorical.data.VGAM) 342 mustart <- mustart[tt.index, 2] 343 y <- y[, 2] 344 } else { 345 } 346 347 348 x.save <- x # Save the original 349 y.save <- y # Save the original 350 w.save <- w # Save the original 351 352 new.coeffs <- .coefstart # Needed for iter = 1 of @weight 353 new.coeffs <- if (length(new.coeffs)) 354 rep_len(new.coeffs, p.lm + plag) else 355 c(rnorm(p.lm, sd = 0.1), rep_len(0, plag)) 356 357 if (!length(etastart)) { 358 etastart <- x[-indices, , drop = FALSE] %*% new.coeffs[1:p.lm] 359 } 360 361 x <- cbind(x, matrix(NA_real_, n, plag)) # Right size now 362 dx <- dimnames(x.save) 363 morenames <- paste("(lag", 1:plag, ")", sep = "") 364 dimnames(x) <- list(dx[[1]], c(dx[[2]], morenames)) 365 366 x <- x[-indices, , drop = FALSE] 367 class(x) <- "matrix" 368 y <- y[-indices] 369 w <- w[-indices] 370 n.save <- n <- n - plag 371 372 more <- vector("list", plag) 373 names(more) <- morenames 374 for (ii in 1:plag) 375 more[[ii]] <- ii + max(unlist(attr(x.save, "assign"))) 376 attr(x, "assign") <- c(attr(x.save, "assign"), more) 377 }), list( .link = link, .p.ar.lag = p.ar.lag, 378 .coefstart = coefstart, .earg = earg ))), 379 linkinv = eval(substitute(function(eta, extra = NULL) { 380 eta2theta(eta, link = .link , earg = .earg) 381 }, list( .link = link, .earg = earg ))), 382 last = eval(substitute(expression({ 383 misc$link <- c(mu = .link ) 384 385 misc$earg <- list(mu = .earg ) 386 387 misc$plag <- plag 388 }), list( .link = link, .earg = earg ))), 389 390 loglikelihood = eval(substitute( 391 function(mu, y, w, residuals = FALSE, eta, extra = NULL, 392 summation = TRUE) { 393 if (residuals) { 394 switch( .link , 395 identitylink = y - mu, 396 loglink = w * (y / mu - 1), 397 reciprocallink = w * (y / mu - 1), 398 inverse = w * (y / mu - 1), 399 w * (y / mu - (1-y) / (1 - mu))) 400 } else { 401 ll.elts <- 402 switch( .link , 403 identitylink = c(w) * (y - mu)^2, 404 loglink = c(w) * (-mu + y * log(mu)), 405 reciprocallink = c(w) * (-mu + y * log(mu)), 406 inverse = c(w) * (-mu + y * log(mu)), 407 c(w) * (y * log(mu) + (1-y) * log1p(-mu))) 408 if (summation) { 409 sum(ll.elts) 410 } else { 411 ll.elts 412 } 413 } 414 }, list( .link = link, .earg = earg ))), 415 416 middle2 = eval(substitute(expression({ 417 realfv <- fv 418 for (ii in 1:plag) { 419 realfv <- realfv + old.coeffs[ii + p.lm] * 420 (x.save[tt.index-ii, 1:p.lm, drop = FALSE] %*% 421 new.coeffs[1:p.lm]) # + 422 } 423 424 true.eta <- realfv + offset 425 mu <- family@linkinv(true.eta, extra) # overwrite mu with correct one 426 }), list( .link = link, .earg = earg ))), 427 vfamily = c("garma", "vglmgam"), 428 validparams = eval(substitute(function(eta, y, extra = NULL) { 429 mu <- eta2theta(eta, link = .link , earg = .earg ) 430 okay1 <- all(is.finite(mu)) 431 okay1 432 }, list( .link = link, .earg = earg ))), 433 deriv = eval(substitute(expression({ 434 dl.dmu <- switch( .link , 435 identitylink = y-mu, 436 loglink = (y - mu) / mu, 437 reciprocallink = (y - mu) / mu, 438 inverse = (y - mu) / mu, 439 (y - mu) / (mu * (1 - mu))) 440 dmu.deta <- dtheta.deta(mu, .link , earg = .earg) 441 Step <- .step # This is another method of adjusting step lengths 442 Step * c(w) * dl.dmu * dmu.deta 443 }), list( .link = link, 444 .step = step, 445 .earg = earg ))), 446 447 weight = eval(substitute(expression({ 448 x[, 1:p.lm] <- x.save[tt.index, 1:p.lm] # Reinstate 449 450 for (ii in 1:plag) { 451 temp <- theta2eta(y.save[tt.index-ii], .link , earg = .earg ) 452 453 454 x[, 1:p.lm] <- x[, 1:p.lm] - 455 x.save[tt.index-ii, 1:p.lm] * new.coeffs[ii + p.lm] 456 x[, p.lm+ii] <- temp - 457 x.save[tt.index-ii, 1:p.lm, drop = FALSE] %*% 458 new.coeffs[1:p.lm] 459 } 460 class(x) <- "matrix" # Added 20020227; 20040226 461 462 if (iter == 1) 463 old.coeffs <- new.coeffs 464 465 X.vlm.save <- lm2vlm.model.matrix(x, Hlist, xij = control$xij) 466 467 vary <- switch( .link , 468 identitylink = 1, 469 loglink = mu, 470 reciprocallink = mu^2, 471 inverse = mu^2, 472 mu * (1 - mu)) 473 c(w) * dtheta.deta(mu, link = .link , earg = .earg )^2 / vary 474 }), list( .link = link, 475 .earg = earg )))) 476} 477 478 479 480 481 482 483 if (FALSE) { 484setClass(Class = "Coef.rrar", representation( 485 "plag" = "integer", 486 "Ranks" = "integer", 487 "omega" = "integer", 488 "C" = "matrix", 489 "D" = "matrix", 490 "H" = "matrix", 491 "Z" = "matrix", 492 "Phi" = "list", # list of matrices 493 "Ak1" = "matrix")) 494 495 496 497Coef.rrar <- function(object, ...) { 498 result = new(Class = "Coef.rrar", 499 "plag" = object@misc$plag, 500 "Ranks" = object@misc$Ranks, 501 "omega" = object@misc$omega, 502 "C" = object@misc$C, 503 "D" = object@misc$D, 504 "H" = object@misc$H, 505 "Z" = object@misc$Z, 506 "Phi" = object@misc$Phi, 507 "Ak1" = object@misc$Ak1) 508} 509 510 511 512 513 514show.Coef.rrar <- function(object) { 515 cat(object@plag) 516} 517 518 519setMethod("Coef", "rrar", 520 function(object, ...) 521 Coef(object, ...)) 522 523 524 525 526setMethod("show", "Coef.rrar", 527 function(object) 528 show.Coef.rrar(object)) 529 530} 531 532 533 534 535 536 537 538 539 540 541dAR1 <- function(x, 542 drift = 0, # Stationarity is the default 543 var.error = 1, ARcoef1 = 0.0, 544 type.likelihood = c("exact", "conditional"), 545 log = FALSE) { 546 547 type.likelihood <- match.arg(type.likelihood, 548 c("exact", "conditional"))[1] 549 550 is.vector.x <- is.vector(x) 551 552 x <- as.matrix(x) 553 drift <- as.matrix(drift) 554 var.error <- as.matrix(var.error) 555 ARcoef1 <- as.matrix(ARcoef1) 556 LLL <- max(nrow(x), nrow(drift), nrow(var.error), nrow(ARcoef1)) 557 UUU <- max(ncol(x), ncol(drift), ncol(var.error), ncol(ARcoef1)) 558 x <- matrix(x, LLL, UUU) 559 drift <- matrix(drift, LLL, UUU) 560 var.error <- matrix(var.error, LLL, UUU) 561 rho <- matrix(ARcoef1, LLL, UUU) 562 563 if (any(abs(rho) > 1)) 564 warning("Values of argument 'ARcoef1' are greater ", 565 "than 1 in absolute value") 566 567 if (!is.logical(log.arg <- log) || length(log) != 1) 568 stop("Bad input for argument 'log'") 569 rm(log) 570 571 ans <- matrix(0.0, LLL, UUU) 572 573 var.noise <- var.error / (1 - rho^2) 574 575 ans[ 1, ] <- dnorm(x = x[1, ], 576 mean = drift[ 1, ] / (1 - rho[1, ]), 577 sd = sqrt(var.noise[1, ]), log = log.arg) 578 ans[-1, ] <- dnorm(x = x[-1, ], 579 mean = drift[-1, ] + rho[-1, ] * x[-nrow(x), ], 580 sd = sqrt(var.error[-1, ]), log = log.arg) 581 582 if (type.likelihood == "conditional") 583 ans[1, ] <- NA 584 585 if (is.vector.x) as.vector(ans) else ans 586} 587 588 589 590if (FALSE) 591AR1.control <- function(epsilon = 1e-6, 592 maxit = 30, 593 stepsize = 1,...){ 594 list(epsilon = epsilon, 595 maxit = maxit, 596 stepsize = stepsize, 597 ...) 598} 599 600 601 602 AR1 <- 603 function(ldrift = "identitylink", 604 lsd = "loglink", 605 lvar = "loglink", 606 lrho = "rhobitlink", 607 idrift = NULL, 608 isd = NULL, 609 ivar = NULL, 610 irho = NULL, 611 imethod = 1, 612 ishrinkage = 0.95, # 0.90; unity means a constant 613 type.likelihood = c("exact", "conditional"), 614 type.EIM = c("exact", "approximate"), 615 var.arg = FALSE, # TRUE, 616 nodrift = FALSE, # TRUE, 617 print.EIM = FALSE, 618 zero = c(if (var.arg) "var" else "sd", "rho") # "ARcoeff1" 619 ) { 620 621 type.likelihood <- match.arg(type.likelihood, 622 c("exact", "conditional"))[1] 623 624 if (length(isd) && !is.Numeric(isd, positive = TRUE)) 625 stop("Bad input for argument 'isd'") 626 627 if (length(ivar) && !is.Numeric(ivar, positive = TRUE)) 628 stop("Bad input for argument 'ivar'") 629 630 if (length(irho) && 631 (!is.Numeric(irho) || any(abs(irho) > 1.0))) 632 stop("Bad input for argument 'irho'") 633 634 type.EIM <- match.arg(type.EIM, c("exact", "approximate"))[1] 635 poratM <- (type.EIM == "exact") 636 637 if (!is.logical(nodrift) || 638 length(nodrift) != 1) 639 stop("argument 'nodrift' must be a single logical") 640 641 if (!is.logical(var.arg) || 642 length(var.arg) != 1) 643 stop("argument 'var.arg' must be a single logical") 644 645 if (!is.logical(print.EIM)) 646 stop("Invalid 'print.EIM'.") 647 648 ismn <- idrift 649 lsmn <- as.list(substitute(ldrift)) 650 esmn <- link2list(lsmn) 651 lsmn <- attr(esmn, "function.name") 652 653 lsdv <- as.list(substitute(lsd)) 654 esdv <- link2list(lsdv) 655 lsdv <- attr(esdv, "function.name") 656 657 lvar <- as.list(substitute(lvar)) 658 evar <- link2list(lvar) 659 lvar <- attr(evar, "function.name") 660 661 lrho <- as.list(substitute(lrho)) 662 erho <- link2list(lrho) 663 lrho <- attr(erho, "function.name") 664 665 n.sc <- if (var.arg) "var" else "sd" 666 l.sc <- if (var.arg) lvar else lsdv 667 e.sc <- if (var.arg) evar else esdv 668 669 new("vglmff", 670 blurb = c(ifelse(nodrift, "Two", "Three"), 671 "-parameter autoregressive process of order-1\n\n", 672 "Links: ", 673 if (nodrift) "" else 674 paste(namesof("drift", lsmn, earg = esmn), ", ", 675 sep = ""), 676 namesof(n.sc , l.sc, earg = e.sc), ", ", 677 namesof("rho", lrho, earg = erho), "\n", 678 "Model: Y_t = drift + rho * Y_{t-1} + error_{t},", 679 "\n", 680 " where 'error_{2:n}' ~ N(0, sigma^2) ", 681 "independently", 682 if (nodrift) ", and drift = 0" else "", 683 "\n", 684 "Mean: drift / (1 - rho)", "\n", 685 "Correlation: rho = ARcoef1", "\n", 686 "Variance: sd^2 / (1 - rho^2)"), 687 688 constraints = eval(substitute(expression({ 689 690 M1 <- 3 - .nodrift 691 dotzero <- .zero 692 # eval(negzero.expression.VGAM) 693 constraints <- 694 cm.zero.VGAM(constraints, x = x, zero = .zero , M = M, 695 predictors.names = predictors.names, 696 M1 = M1) 697 698 }), list( .zero = zero, 699 .nodrift = nodrift ))), 700 701 infos = eval(substitute(function(...) { 702 list(M1 = 3 - .nodrift , 703 Q1 = 1, 704 expected = TRUE, 705 multipleResponse = TRUE, 706 type.likelihood = .type.likelihood , 707 ldrift = if ( .nodrift ) NULL else .lsmn , 708 edrift = if ( .nodrift ) NULL else .esmn , 709 lvar = .lvar , 710 lsd = .lsdv , 711 evar = .evar , 712 esd = .esdv , 713 lrho = .lrho , 714 erho = .erho , 715 zero = .zero ) 716 }, list( .lsmn = lsmn, .lvar = lvar, .lsdv = lsdv, .lrho = lrho, 717 .esmn = esmn, .evar = evar, .esdv = esdv, .erho = erho, 718 .type.likelihood = type.likelihood, 719 .nodrift = nodrift, .zero = zero))), 720 721 initialize = eval(substitute(expression({ 722 extra$M1 <- M1 <- 3 - .nodrift 723 check <- w.y.check(w = w, y = y, 724 Is.positive.y = FALSE, 725 ncol.w.max = Inf, 726 ncol.y.max = Inf, 727 out.wy = TRUE, 728 colsyperw = 1, 729 maximize = TRUE) 730 w <- check$w 731 y <- check$y 732 if ( .type.likelihood == "conditional") { 733 w[1, ] <- 1.0e-6 734 } else { 735 if (!(.nodrift )) 736 w[1, ] <- 1.0e-1 737 } 738 739 NOS <- ncoly <- ncol(y) 740 n <- nrow(y) 741 M <- M1*NOS 742 743 var.names <- param.names("var", NOS, skip1 = TRUE) 744 sdv.names <- param.names("sd", NOS, skip1 = TRUE) 745 smn.names <- if ( .nodrift ) NULL else 746 param.names("drift", NOS, skip1 = TRUE) 747 rho.names <- param.names("rho", NOS, skip1 = TRUE) 748 749 mynames1 <- smn.names 750 mynames2 <- if ( .var.arg ) var.names else sdv.names 751 mynames3 <- rho.names 752 753 predictors.names <- 754 c(if ( .nodrift ) NULL else 755 namesof(smn.names, .lsmn , earg = .esmn , tag = FALSE), 756 if ( .var.arg ) 757 namesof(var.names, .lvar , earg = .evar , tag = FALSE) else 758 namesof(sdv.names, .lsdv , earg = .esdv , tag = FALSE), 759 namesof(rho.names, .lrho , earg = .erho , tag = FALSE)) 760 761 predictors.names <- predictors.names[interleave.VGAM(M, M1 = M1)] 762 763 if ( .nodrift ) 764 y <- scale(y, scale = FALSE) 765 766 if (!length(etastart)) { 767 init.smn <- Init.mu(y = y, w = w, imethod = .imethod , # x = x, 768 imu = .ismn , ishrinkage = .ishrinkage , 769 pos.only = FALSE) 770 771 init.rho <- matrix(if (length( .irho )) .irho else 0.1, 772 n, NOS, byrow = TRUE) 773 init.sdv <- matrix(if (length( .isdv )) .isdv else 1.0, 774 n, NOS, byrow = TRUE) 775 init.var <- matrix(if (length( .ivar )) .ivar else 1.0, 776 n, NOS, byrow = TRUE) 777 for (jay in 1: NOS) { 778 mycor <- cor(y[-1, jay], y[-n, jay]) 779 init.smn[ , jay] <- mean(y[, jay]) * (1 - mycor) 780 if (!length( .irho )) 781 init.rho[, jay] <- sign(mycor) * min(0.95, abs(mycor)) 782 if (!length( .ivar )) 783 init.var[, jay] <- var(y[, jay]) * (1 - mycor^2) 784 if (!length( .isdv )) 785 init.sdv[, jay] <- sqrt(init.var[, jay]) 786 } # for 787 788 etastart <- 789 cbind(if ( .nodrift ) NULL else 790 theta2eta(init.smn, .lsmn , earg = .esmn ), 791 if ( .var.arg ) 792 theta2eta(init.var, .lvar , earg = .evar ) else 793 theta2eta(init.sdv, .lsdv , earg = .esdv ), 794 theta2eta(init.rho, .lrho , earg = .erho )) 795 796 etastart <- etastart[, interleave.VGAM(M, M1 = M1), drop = FALSE] 797 } # end of etastart 798 799 }), list( .lsmn = lsmn, .lrho = lrho, .lsdv = lsdv, .lvar = lvar, 800 .esmn = esmn, .erho = erho, .esdv = esdv, .evar = evar, 801 .ismn = ismn, .irho = irho, .isdv = isd , .ivar = ivar, 802 .type.likelihood = type.likelihood, 803 .ishrinkage = ishrinkage, .poratM = poratM, 804 .var.arg = var.arg, 805 .nodrift = nodrift, 806 .imethod = imethod ))), 807 808 linkinv = eval(substitute(function(eta, extra = NULL) { 809 M1 <- 3 - .nodrift 810 NOS <- ncol(eta)/M1 811 ar.smn <- if ( .nodrift ) 0 else 812 eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 813 .lsmn , earg = .esmn ) 814 ar.rho <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 815 .lrho , earg = .erho ) 816 ar.smn / (1 - ar.rho) 817 818 }, list ( .lsmn = lsmn, .lrho = lrho , .lsdv = lsdv, .lvar = lvar , 819 .var.arg = var.arg, .type.likelihood = type.likelihood, 820 .nodrift = nodrift, 821 .esmn = esmn, .erho = erho , .esdv = esdv, .evar = evar ))), 822 823 last = eval(substitute(expression({ 824 if (any(abs(ar.rho) > 1)) 825 warning("Regularity conditions are violated at the final", 826 "IRLS iteration, since 'abs(rho) > 1") 827 828 M1 <- extra$M1 829 830 temp.names <- c(mynames1, mynames2, mynames3) 831 temp.names <- temp.names[interleave.VGAM(M1 * ncoly, M1 = M1)] 832 833 misc$link <- rep_len( .lrho , M1 * ncoly) 834 misc$earg <- vector("list", M1 * ncoly) 835 names(misc$link) <- 836 names(misc$earg) <- temp.names 837 for (ii in 1:ncoly) { 838 if ( !( .nodrift )) 839 misc$link[ M1*ii-2 ] <- .lsmn 840 misc$link[ M1*ii-1 ] <- if ( .var.arg ) .lvar else .lsdv 841 misc$link[ M1*ii ] <- .lrho 842 if ( !( .nodrift )) 843 misc$earg[[M1*ii-2]] <- .esmn 844 misc$earg[[M1*ii-1]] <- if ( .var.arg ) .evar else .esdv 845 misc$earg[[M1*ii ]] <- .erho 846 } 847 848 misc$type.likelihood <- .type.likelihood 849 misc$var.arg <- .var.arg 850 misc$M1 <- M1 851 misc$expected <- TRUE 852 misc$imethod <- .imethod 853 misc$multipleResponses <- TRUE 854 misc$nodrift <- .nodrift 855 856 }), list( .lsmn = lsmn, .lrho = lrho, .lsdv = lsdv, .lvar = lvar, 857 .esmn = esmn, .erho = erho, .esdv = esdv, .evar = evar, 858 .irho = irho, .isdv = isd , .ivar = ivar, 859 .nodrift = nodrift, .poratM = poratM, 860 .var.arg = var.arg, .type.likelihood = type.likelihood, 861 .imethod = imethod ))), 862 863 loglikelihood = eval(substitute( 864 function(mu, y, w, residuals= FALSE, eta, 865 extra = NULL, summation = TRUE) { 866 M1 <- 3 - .nodrift 867 NOS <- ncol(eta)/M1 868 869 if ( .var.arg ) { 870 ar.var <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 871 .lvar , earg = .evar ) 872 ar.sdv <- sqrt(ar.var) 873 } else { 874 ar.sdv <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 875 .lsdv , earg = .esdv ) 876 ar.var <- ar.sdv^2 877 } 878 ar.smn <- if ( .nodrift ) 0 else 879 eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 880 .lsmn , earg = .esmn ) 881 ar.rho <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 882 .lrho , earg = .erho ) 883 884 if (residuals) { 885 stop("Loglikelihood not implemented yet to handle", 886 "residuals.") 887 } else { 888 loglik.terms <- 889 c(w) * dAR1(x = y, 890 drift = ar.smn, 891 var.error = ar.var, 892 type.likelihood = .type.likelihood , 893 ARcoef1 = ar.rho, log = TRUE) 894 loglik.terms <- as.matrix(loglik.terms) 895 896 if (summation) { 897 sum(if ( .type.likelihood == "exact") loglik.terms else 898 loglik.terms[-1, ] ) 899 } else { 900 loglik.terms 901 } 902 } 903 904 }, list( .lsmn = lsmn, .lrho = lrho , .lsdv = lsdv, .lvar = lvar , 905 .var.arg = var.arg, .type.likelihood = type.likelihood, 906 .nodrift = nodrift, 907 .esmn = esmn, .erho = erho , 908 .esdv = esdv, .evar = evar ))), 909 910 vfamily = c("AR1"), 911 validparams = eval(substitute(function(eta, y, extra = NULL) { 912 M1 <- 3 - .nodrift 913 n <- nrow(eta) 914 NOS <- ncol(eta)/M1 915 ncoly <- NCOL(y) 916 917 if ( .var.arg ) { 918 ar.var <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 919 .lvar , earg = .evar ) 920 ar.sdv <- sqrt(ar.var) 921 } else { 922 ar.sdv <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 923 .lsdv , earg = .esdv ) 924 ar.var <- ar.sdv^2 925 } 926 ar.smn <- if ( .nodrift ) matrix(0, n, NOS) else 927 eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 928 .lsmn , earg = .esmn ) 929 ar.rho <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 930 .lrho , earg = .erho ) 931 okay1 <- all(is.finite(ar.sdv)) && all(0 < ar.sdv) && 932 all(is.finite(ar.smn)) && 933 all(is.finite(ar.rho)) 934 okay1 935 }, list( .lsmn = lsmn, .lrho = lrho , .lsdv = lsdv, .lvar = lvar , 936 .var.arg = var.arg, .type.likelihood = type.likelihood, 937 .nodrift = nodrift, 938 .esmn = esmn, .erho = erho , 939 .esdv = esdv, .evar = evar ))), 940 941 simslot = eval(substitute( 942 function(object, nsim) { 943 944 pwts <- if (length(pwts <- object@prior.weights) > 0) 945 pwts else weights(object, type = "prior") 946 if (any(pwts != 1)) 947 warning("ignoring prior weights") 948 eta <- predict(object) 949 fva <- fitted(object) 950 M1 <- 3 - .nodrift 951 NOS <- ncol(eta)/M1 952 953 if ( .var.arg ) { 954 ar.var <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 955 .lvar , earg = .evar ) 956 ar.sdv <- sqrt(ar.var) 957 } else { 958 ar.sdv <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 959 .lsdv , earg = .esdv ) 960 ar.var <- ar.sdv^2 961 } 962 ar.smn <- if ( .nodrift ) matrix(0, n, NOS) else 963 eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 964 .lsmn , earg = .esmn ) 965 ar.rho <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 966 .lrho , earg = .erho ) 967 968 ans <- array(0, c(nrow(eta), NOS, nsim)) 969 for (jay in 1:NOS) { 970 ans[1, jay, ] <- rnorm(nsim, m = fva[1, jay], # zz 971 sd = sqrt(ar.var[1, jay])) 972 for (ii in 2:nrow(eta)) 973 ans[ii, jay, ] <- ar.smn[ii, jay] + 974 ar.rho[ii, jay] * ans[ii-1, jay, ] + 975 rnorm(nsim, sd = sqrt(ar.var[ii, jay])) 976 } 977 ans <- matrix(c(ans), c(nrow(eta) * NOS, nsim)) 978 ans 979 980 }, list( .lsmn = lsmn, .lrho = lrho , .lsdv = lsdv, .lvar = lvar , 981 .var.arg = var.arg, .type.likelihood = type.likelihood, 982 .nodrift = nodrift, 983 .esmn = esmn, .erho = erho , 984 .esdv = esdv, .evar = evar ))), 985 986 987 deriv = eval(substitute(expression({ 988 M1 <- 3 - .nodrift 989 NOS <- ncol(eta)/M1 990 ncoly <- NCOL(y) 991 992 if ( .var.arg ) { 993 ar.var <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 994 .lvar , earg = .evar ) 995 ar.sdv <- sqrt(ar.var) 996 } else { 997 ar.sdv <- eta2theta(eta[, M1*(1:NOS) - 1, drop = FALSE], 998 .lsdv , earg = .esdv ) 999 ar.var <- ar.sdv^2 1000 } 1001 1002 ar.smn <- if ( .nodrift ) matrix(0, n, NOS) else 1003 eta2theta(eta[, M1*(1:NOS) - 2, drop = FALSE], 1004 .lsmn , earg = .esmn ) 1005 1006 ar.rho <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE], 1007 .lrho , earg = .erho ) 1008 1009 if (any(abs(ar.rho) < 1e-2)) 1010 warning("Estimated values of 'rho' are too close to zero.") 1011 1012 help2 <- (length(colnames(x)) >= 2) 1013 myMeans <- matrix(colMeans(y), nrow = n, ncol = NOS, by = TRUE) 1014 yLag <- matrix(y, ncol = NOS) 1015 temp4 <- matrix(0.0, nrow = n, ncol = NOS) 1016 temp4[-1, ] <- y[-1, , drop = FALSE] - ar.smn[-1, , drop = FALSE] 1017 yLag[-1, ] <- y[-n, ] 1018 1019 temp1 <- matrix(0.0, nrow = n, ncol = NOS) 1020 temp1[-1, ] <- y[-1, , drop = FALSE] - (ar.smn[-1, ,drop = FALSE] + 1021 ar.rho[-1, , drop = FALSE] * y[-n, , drop = FALSE]) 1022 temp1[1, ] <- y[1, ] - ar.smn[1, ] 1023 dl.dsmn <- temp1 / ar.var 1024 dl.dsmn[1, ] <- ( (y[1, ] - myMeans[1, ]) * 1025 (1 + ar.rho[1, ]) ) / ar.var[1, ] 1026 1027 if ( .var.arg ) { 1028 dl.dvarSD <- temp1^2 / ( 2 * ar.var^2) - 1 / (2 * ar.var) 1029 dl.dvarSD[1, ] <- ( (1 - ar.rho[1, ]^2) * (y[1, ] - 1030 myMeans[1, ])^2 ) /(2 * ar.var[1, ]^2) - 1 / (2 * ar.var[1, ]) 1031 } else { 1032 dl.dvarSD <- temp1^2 / ar.sdv^3 - 1 / ar.sdv 1033 dl.dvarSD[1, ] <- ( (1 - ar.rho[1, ]^2) * 1034 (y[1, ] - myMeans[1, ])^2 ) / ar.sdv[1, ]^3 - 1/ar.sdv[1, ] 1035 } 1036 1037 dl.drho <- rbind(rep_len(0, 1), 1038 ( (y[-n, , drop = FALSE] - myMeans[-n, ]) * 1039 temp1[-1, , drop = FALSE ] )/ ar.var[-1, ] ) 1040 dl.drho[1, ] <- (ar.rho[1, ] * (y[1, ] - myMeans[1, ])^2 ) / 1041 ar.var[1, ] - ar.rho[1, ] / (1 - ar.rho[1, ]^2) 1042 1043 dsmn.deta <- dtheta.deta(ar.smn, .lsmn , earg = .esmn ) 1044 drho.deta <- dtheta.deta(ar.rho, .lrho , earg = .erho ) 1045 if ( .var.arg ) { 1046 dvarSD.deta <- dtheta.deta(ar.var, .lvar , earg = .evar ) 1047 } else { 1048 dvarSD.deta <- dtheta.deta(ar.sdv, .lsdv , earg = .esdv ) 1049 } 1050 1051 myderiv <- 1052 c(w) * cbind(if ( .nodrift ) NULL else dl.dsmn * dsmn.deta, 1053 dl.dvarSD * dvarSD.deta, 1054 dl.drho * drho.deta) 1055 myderiv <- myderiv[, interleave.VGAM(M, M1 = M1)] 1056 myderiv 1057 1058 }), list( .lsmn = lsmn, .lrho = lrho, .lsdv = lsdv, .lvar = lvar, 1059 .esmn = esmn, .erho = erho, .esdv = esdv, .evar = evar, 1060 .nodrift = nodrift , 1061 .var.arg = var.arg, 1062 .type.likelihood = type.likelihood ))), 1063 1064 weight = eval(substitute(expression({ 1065 1066 ned2l.dsmn <- 1 / ar.var 1067 ned2l.dsmn[1, ] <- ( (1 + ar.rho[1, ]) / (1 - ar.rho[1, ]) ) * 1068 (1 / ar.var[1, ]) 1069 # Here, same results for the first and t > 1 observations. 1070 ned2l.dvarSD <- if ( .var.arg ) 1 / (2 * ar.var^2) else 2 / ar.var 1071 gamma0 <- (1 - help2) * ar.var/(1 - ar.rho^2) + 1072 help2 * (yLag - myMeans)^2 1073 ned2l.drho <- gamma0 / ar.var 1074 ned2l.drho[1, ] <- 2 * ar.rho[1, ]^2 / (1 - ar.rho[1, ]^2)^2 1075 ned2l.drdv <- matrix(0.0, nrow = n, ncol = NOS) 1076 ned2l.drdv[1, ] <- 2 * temp4[1, ] / 1077 ((1 - temp4[1, ]^2) * ar.sdv[1, ]) 1078 ncol.wz <- M + (M - 1) + ifelse( .nodrift , 0, M - 2) 1079 ncol.pf <- 3 * (M + ( .nodrift ) ) - 3 1080 wz <- matrix(0, nrow = n, ncol = ncol.wz) 1081 helpPor <- .poratM 1082 1083 pf.mat <- if (helpPor) 1084 AR1EIM(x = scale(y, scale = FALSE), 1085 var.arg = .var.arg , 1086 p.drift = 0, 1087 WNsd = ar.sdv, 1088 ARcoeff1 = ar.rho ) else 1089 array(0.0, dim= c(n, NOS, ncol.pf)) 1090 1091 if (!( .nodrift )) 1092 wz[, M1*(1:NOS) - 2] <- ( (helpPor) * pf.mat[, , 1] + 1093 (1 - (helpPor)) * ned2l.dsmn) * dsmn.deta^2 1094 wz[, M1*(1:NOS) - 1] <- ( (helpPor) * pf.mat[, , 2 ] + 1095 (1 - (helpPor)) * ned2l.dvarSD) * dvarSD.deta^2 1096 wz[, M1*(1:NOS) ] <- ( (helpPor) * pf.mat[, , 3] + 1097 (1 - (helpPor)) * ned2l.drho) * drho.deta^2 1098 wz[, M1*(1:NOS) + (M - 1) ] <- ((helpPor) * pf.mat[, , 4] + 1099 (1 - (helpPor)) * ned2l.drdv) * drho.deta * dvarSD.deta 1100 1101 wz <- w.wz.merge(w = w, wz = wz, n = n, 1102 M = ncol.wz, ndepy = NOS) 1103 1104 if ( .print.EIM ) { 1105 wz2 <- matrix(0, nrow = n, ncol = ncol.wz) 1106 if (!(.nodrift )) 1107 wz2[, M1*(1:NOS) - 2] <- ned2l.dsmn 1108 wz2[, M1*(1:NOS) - 1] <- 1109 if ( .var.arg ) 1 / (2 * ar.var^2) else 2 / ar.var 1110 wz2[, M1*(1:NOS) ] <- ned2l.drho 1111 1112 wz2 <- wz2[, interleave.VGAM( M1 * NOS, M1)] 1113 if (NOS > 1) { 1114 1115 matAux1 <- matAux2 <- matrix(NA_real_, nrow = n, ncol = NOS) 1116 approxMat <- array(wz2[, 1:(M1*NOS)], dim = c(n, M1, NOS)) 1117 for (kk in 1:NOS) { 1118 matAux1[, kk] <- rowSums(approxMat[, , kk]) 1119 matAux2[, kk] <- rowSums(pf.mat[, kk , ]) 1120 } 1121 matAux <- cbind(matAux1, if (.poratM ) matAux2 else NULL) 1122 colnames(matAux) <- c(param.names("ApproxEIM.R", NOS), 1123 if (!(.poratM )) NULL else 1124 param.names("ExactEIM.R", NOS)) 1125 1126 matAux <- matAux[, interleave.VGAM( (1 + .poratM) * NOS, 1127 M1 = 1 + .poratM)] 1128 } else { 1129 1130 matAux <- cbind(rowSums(wz2), 1131 if (helpPor) 1132 rowSums(pf.mat[, 1, ][, 1:3]) else NULL) 1133 colnames(matAux) <- c("Approximate", 1134 if (helpPor) "Exact" else NULL) 1135 1136 } 1137 print(matAux[1:10, , drop = FALSE]) 1138 } 1139 1140 wz 1141 1142 }), list( .var.arg = var.arg, .type.likelihood = type.likelihood, 1143 .nodrift = nodrift, .poratM = poratM, 1144 .print.EIM = print.EIM ))) 1145 ) 1146 } 1147 1148 1149 1150 1151 1152 1153