#' AR steady state mean and variance #' #' Finds mean and variance after many periods, steady state only if roots stable. #' #' @param A Either a matrix coefficient on a single lag, or an array of matrix #' coefficients on several lags. #' @param sigma covariance matrix of system disturbances #' @param const vector of constant terms, or NULL if none #' @param y0 initial values for y in the iteration. Usually just a zero vector #' @param log2T Log to base 2 of the number of iterations of the doubling algorithm #' Increase above the default value if the system is in fact stable and/or #' you really want to be closer to steady state or to depress likelihood #' for a non-stationary model. #' @details Computes the conditional distribution of \code{y, 2^{log2T}} periods #' after an artificial initial date. #' @return a list with elements #' \describe{ #' \item{\code{Omega}}{the covariance matrix} #' \item{\code{mu}{the mean} #' \item{\code{B}{A matrix, all of whose elements (except the last column if there is a constant) #' should be very small if the result is in fact close to the steady state. If not, #' you might want to increase log2T} #' } #' @export doubleic <- function(A, sigma, const=c(0, NULL), y0=c(0,0), log2T=8) { ## if there is no constant, it would be more efficient to use a routine that did not ## compute the mean along with the covariance matrix. ## To start from 0 initial condition, but with a constant, set y0 to a vector of zeros. ## To invoke recursively with output from previous run, const=NULL, y0=mu from previous ## returned value, if ( is.null(dim(A))) { if (length(A) == 1) { A <- array(A, c(1,1,1)) } else { A <- array(A, c(1, 1, length(A))) } } if (length(dim(A)) == 2) #first order system A <- array(A, c(dim(A), 1)) ## assumes A is now an n by n by nlags array ny <- dim(A)[1] nnl <- prod(dim(A)[2:3]) if (!is.null(const) ) { work <- matrix(0, nnl + 1, nnl + 1) } else { work <- matrix(0, nnl, nnl) } work[1:ny, 1:ny] <- sigma sigma <- work A <- sysmat(A, const) if (!is.null(const)) mu <- c(y0, 1) else mu <- y0 B <- A Omega <- sigma for (nit in 1:log2T) { Omega <- Omega + B %*% Omega %*% t(B) B <- B %*% B } mu <- B %*% mu ## If one does not want to commit to a log2T value, this routine can be called recursively, ## with A, sigma const arguments replaced by the returned B, Omega, mu, and the calculation ## stopped when norm of upper left block (excluding constants) of $B$ becomes small enough. return(list(Omega=Omega, mu=mu, B=B )) }