estVARMA <- function(ydata, By0, Be0, Bx0, Sig0, H0=NULL, useIC=FALSE, nit) { ny <- dim(By0)[1] nly <- dim(By0)[3] nle <- dim(Be0)[3] ne <- dim(Be0)[2] nby <- ny * ny * nly nbe <- ny * ne * nle nbx <- ny nsig <- ny * (ny + 1)/2 ## lines below condition on initial y's and assume initial eps is 0. ## An alternative would be to use steady state distribution, use all data in lh. nstate <- ny * nly + ne * nle + 1 wrapper <- function(param) { By <- array(param[1:nby], c(ny,ny,nly)) Be <- array(param[nby + 1:nbe], c(ny, ne, nle)) Bx <- param[nby + nbe + 1:nbx] Sig <- matrix(0, ny, ny) Sig[upper.tri(Sig, diag=TRUE)] <- param[nby + nbe + nbx + 1:nsig] Sig[lower.tri(Sig)] <- t(Sig)[lower.tri(Sig)] mats <- armasysmat(By, Be, Bx, Sig) G <- mats$sysmat M <- cschol(mats$sigkf) H <- diag(1, dim(By)[1], dim(G)[2]) #observed y's at top of state. if (useIC) { icmv <- doubleic(G, mats$sigkf, const=NULL, y0=c(rep(0, dim(G)[1]-1), 1)) icm <- icmv$mu icsig <- icmv$Omega ##browser() vml <- VARMAlh(ydata, G, H, M, icm, icsig) } else { icm <- c(t(ydata[nly:1, ]), rep(0, ne * nle), 1) icsig <- matrix(0, nstate, nstate) ##browser() vml <- VARMAlh(ydata[-(1:nly), ], G, H, M, icm, icsig) } return(-sum(vml$lh)) } if (is.null(H0)) H0 <- diag(nby + nbe + nbx + nsig) * 1e-6 csout <- csminwelNew(wrapper, c(By0,Be0,Bx0, Sig0[upper.tri(Sig0, diag=TRUE)]), H0, nit=nit) return(csout) }