#######################################################################################################
# Parameter estimation using R's optim function
# Minimal error checking is done. You should run is.marssMLE() before calling this.
# Likelihood computation is via Kalman filter
#######################################################################################################
MARSSoptim = function(MLEobj) {
# This function does not check if user specified a legal MLE object.
tmp = is.marssMLE(MLEobj)
if(!isTRUE(tmp)) {
cat(tmp)
stop("Stopped in MARSSoptim() because marssMLE object is incomplete or inconsistent.\n", call.=FALSE)
}
for(elem in c("Q","R")){
if(dim(MLEobj$model$free[[elem]])[3]>1)
stop(paste("Stopped in MARSSoptim() because function does not allow estimated part of ",elem," to be time-varying.\n",sep=""), call.=FALSE)
}
#the is.marssMODEL code has an additional test that the blocks are diagonal or unconstrained in the varcov matrices
# if(MLEobj$model$diffuse)
# stop(paste("Stopped in MARSSoptim() because version 3.0 does not allow diffuse=TRUE.\n",sep=""), call.=FALSE)
## attach would be risky here since user might have one of these variables in their workspace
modelObj=MLEobj[["marss"]]
y = modelObj$data #must have time going across columns
free = modelObj$free
fixed = modelObj$fixed
tmp.inits = MLEobj$start
control=MLEobj$control
par.dims=attr(modelObj,"model.dims")
m=par.dims[["x"]][1]
n=par.dims[["y"]][1]
## Set up the control list for optim; only pass in optim control elements
control.names=c("trace", "fnscale", "parscale", "ndeps", "maxit", "abstol", "reltol", "alpha", "beta", "gamma", "REPORT", "type", "lmm", "factr","pgtol", "temp", "tmax")
optim.control=list()
for(elem in control.names){
if(!is.null(control[[elem]])) optim.control[[elem]]=control[[elem]]
}
if(is.null(control[["lower"]])){ lower = -Inf
}else lower=control[["lower"]]
if(is.null(control[["upper"]])){ upper = Inf
}else upper=control$upper
if(control$trace==-1)control$trace=0
#The code is used to set things up to use MARSSvectorizeparam to just select inits for the estimated parameters
#Q=t(chol(Q)%*%chol(Q)); t(chol(Q)) has 0 in upper triangle
tmp.MLEobj = MLEobj
#This is needed for the likelihood calculation
tmp.MLEobj$fixed.original = tmp.MLEobj$marss$fixed
tmp.MLEobj$free.original = tmp.MLEobj$marss$free
tmp.MLEobj$par = tmp.inits #set initial conditions for estimated parameters
for(elem in c("Q","R","V0")){ #need the chol for these
d=sub3D(free[[elem]],t=1) #free[[elem]] is required to be time constant
f=sub3D(fixed[[elem]],t=1)
the.par=unvec(f+d%*%tmp.inits[[elem]], dim=par.dims[[elem]][1:2])
is.zero=diag(the.par)==0 #where the 0s on diagonal are
if(any(is.zero)) diag(the.par)[is.zero]=1 #so the chol doesn't fail if there are zeros on the diagonal
the.par=t(chol(the.par)) #convert to transpose of chol
if(any(is.zero)) diag(the.par)[is.zero]=0 #set back to 0
if(!is.fixed(free[[elem]])){
#from f+Dm=M so m = solve(crossprod(d))%*%t(d)%*%(vec(the.par)-f)
#but if d!=0,then f==0. if f!-0, then d==0.
#Thus crossprod(d))%*%t(d) has 0 cols where fs appear in the.par and f is not needed
tmp.MLEobj$par[[elem]] = solve(crossprod(d))%*%t(d)%*%vec(the.par)
}else{ tmp.MLEobj$par[[elem]] = matrix(0,0,1) }
#when being passed to optim, pars for var-cov mat is the chol, so need to reset free so we can get the L=t(chol) matrix
#we don't need to reset fixed because it won't be used;
#step 1, compute the D matrix corresponding to upper.tri=0 in t(chol)
#note this only works because it is required that
#a) if f!=0, d=0. so 1+a never appears in var-cov mat b) a+b never appears in a var-cov mat, c) for BFGS, var-cov mat is time-invariant
tmp.list.mat=fixed.free.to.formula(f,d, par.dims[[elem]][1:2])
tmp.list.mat[upper.tri(tmp.list.mat)]=0 #set upper tri to zero
tmp.MLEobj$marss$free[[elem]]=convert.model.mat(tmp.list.mat)$free
}
# will return the inits only for the estimated parameters
pars = MARSSvectorizeparam(tmp.MLEobj)
if(substr(tmp.MLEobj$method,1,4)=="BFGS"){ optim.method="BFGS" }else{ optim.method="something wrong" }
kf.function=MLEobj$fun.kf #used for printing
optim.output = try(optim(pars, neglogLik, MLEobj=tmp.MLEobj, method = optim.method, lower = lower, upper = upper, control = optim.control, hessian = FALSE), silent=TRUE )
if(class(optim.output)=="try-error"){ #try MARSSkfss if the user did not use it
if( MLEobj$fun.kf!="MARSSkfss" ){ #if user did not request MARSSkf
cat("MARSSkfas returned error. Trying MARSSkfss.\n")
tmp.MLEobj$fun.kf="MARSSkfss"
kf.function="MARSSkfss" #used for printing
optim.output = try(optim(pars, neglogLik, MLEobj=tmp.MLEobj, method = optim.method, lower = lower, upper = upper, control = optim.control, hessian = FALSE), silent=TRUE )
}
}
#error returned
if(class(optim.output)=="try-error"){
optim.output = list(convergence=53, message=c("MARSSkfas and MARSSkfss tried to compute log likelihood and encountered numerical problems.\n", sep=""))
}
MLEobj.return=MLEobj
MLEobj.return$iter.record=optim.output$message
# MLEobj.return$control=MLEobj$control
# MLEobj.return$model=MLEobj$model
MLEobj.return$start = tmp.inits #set to what was used here
MLEobj.return$convergence = optim.output$convergence
if(optim.output$convergence %in% c(1,0)) {
if((!control$silent || control$silent==2) && optim.output$convergence==0) cat(paste("Success! Converged in ",optim.output$counts[1]," iterations.\n","Function ",kf.function," used for likelihood calculation.\n",sep=""))
if((!control$silent || control$silent==2) && optim.output$convergence==1) cat(paste("Warning! Max iterations of ", control$maxit," reached before convergence.\n","Function ", kf.function, " used for likelihood calculation.\n", sep=""))
tmp.MLEobj = MARSSvectorizeparam(tmp.MLEobj, optim.output$par)
#par has the fixed and estimated values using t chol of Q and R
#back transform Q, R and V0 if needed from chol form to usual form
for(elem in c("Q","R","V0")){ #this works because by def fixed and free blocks of var-cov mats are independent
if(!is.fixed(modelObj$free[[elem]])) #get a new par if needed
{
d=sub3D(tmp.MLEobj$marss$free[[elem]],t=1) #this will be the one with the upper tri zero-ed out but ok since symmetric
par.dim=par.dims[[elem]][1:2]
L=unvec(tmp.MLEobj$marss$free[[elem]][,,1]%*%tmp.MLEobj$par[[elem]],dim=par.dim) #this by def will have 0 row/col at the fixed values
the.par = tcrossprod(L)#L%*%t(L)
tmp.MLEobj$par[[elem]]=solve(crossprod(d))%*%t(d)%*%vec(the.par)
}
} #end for
pars = MARSSvectorizeparam(tmp.MLEobj) #now the pars values have been adjusted back to normal scaling
#now put the estimated values back into the original MLEobj; fixed and free matrices as in original
MLEobj.return = MARSSvectorizeparam(MLEobj.return, pars)
kf.out = MARSSkf(MLEobj.return)
}else{
if(optim.output$convergence==10) optim.output$message=c("degeneracy of the Nelder-Mead simplex\n",paste("Function ",kf.function," used for likelihood calculation.\n",sep=""),optim.output$message)
optim.output$counts = NULL
if( !control$silent ) cat("MARSSoptim() stopped with errors. No parameter estimates returned.\n")
if( control$silent==2 ) cat("MARSSoptim() stopped with errors. No parameter estimates returned. See $errors in output for details.\n")
MLEobj.return$par = NULL
MLEobj.return$errors = optim.output$message
kf.out = NULL
}
if(!is.null(kf.out)){
if(control$trace>0) MLEobj.return$kf = kf.out
MLEobj.return$states = kf.out$xtT
MLEobj.return$numIter = optim.output$counts[1]
MLEobj.return$logLik = kf.out$logLik
}
MLEobj.return$method = MLEobj$method
## Add AIC and AICc to the object
if(!is.null(kf.out)) MLEobj.return = MARSSaic(MLEobj.return)
# I don't want these added here
# ## Calculate confidence intervals based on state std errors, see caption of Fig 6.3 (p337) Shumway and Stoffer
# if(!is.null(kf.out)){
# TT = dim(MLEobj.return$marss$data)[2]; m = dim(MLEobj.return$marss$fixed$x0)[1]
# if(m == 1) states.se = sqrt(matrix(kf.out$VtT[,,1:TT], nrow=1))
# if(m > 1) {
# states.se = matrix(0, nrow=m, ncol=TT)
# for(i in 1:TT) states.se[,i] = t(sqrt(takediag(kf.out$VtT[,,i])))
# }
# MLEobj.return$states.se = states.se
# }
return(MLEobj.return)
}
neglogLik = function(x, MLEobj=NULL){ #NULL assignment needed for optim call syntax
#MLEobj is tmp.MLEobj so has altered free and fixed
#x is the paramvector
#update the MLEobj by putting the estimated pars from optim in
MLEobj = MARSSvectorizeparam(MLEobj, x)
free=MLEobj$marss$free
pars=MLEobj$par
par.dims=attr(MLEobj[["marss"]],"model.dims")
for(elem in c("Q","R","V0")){
if(!is.fixed(free[[elem]])) #recompute par if needed since par in parlist is transformed
{
d=sub3D(free[[elem]],t=1) #this will be the one with the upper tri zero-ed out
par.dim=par.dims[[elem]][1:2]
#t=1 since D not allowed to be time-varying; since code 4 lines down won't work otherwise
L=unvec(d%*%pars[[elem]],dim=par.dim) #this by def will have 0 row/col at the fixed values
the.par = tcrossprod(L)#L%*%t(L)
#from f+Dm=M and if f!=0, D==0 so can leave off f
MLEobj$par[[elem]]=solve(crossprod(d))%*%t(d)%*%vec(the.par)
#solve(t(d)%*%d)%*%t(d)%*%vec(the.par)
}
} #end for over elem
#This function is passed a special MLEobj with a marss.original element
MLEobj$marss$fixed = MLEobj$fixed.original
MLEobj$marss$free = MLEobj$free.original
#kfsel selects the Kalman filter / smoother function based on MLEobj$fun.kf
negLL = MARSSkf( MLEobj, only.logLik=TRUE, return.lag.one=FALSE )$logLik
-1*negLL
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.