Unable to correctly estimate ARMA parameters

[ypwong22]

Dear Dr. Helske,
I am having trouble using the KFAS package to estimate ARMA coefficients. I did a test with the Nile river data, and the estimated coefficients differ a lot from what can be estimated by the standard arima in R. The codes are reproduced below:

data(Nile)
Nile = diff(Nile)
model = SSModel( Nile ~ SSMarima(ar=c(0.1,0.2),ma=c(0.1)), H=1e-7 )
fn <- function(pars,model,...){

update T (which contains AR coefficients) and R (which contains MA coefficients)

model$T[2,2,] <- pars[1]
model$T[3,2,] <- pars[2]
model$R[3,,1] <- pars[3]

update the ARIMA Q

model$Q[1,1,1] <- exp( pars[4] )

model
}
inits = c(0.1,0.5,0.5,100)
modelFit<-fitSSM( model=model, inits = inits,
updatefn = fn, method='BFGS', control=list(REPORT=1,trace=1) )

the above does not coincide with the estimation result of

f = arima(x =Nile, order = c(2,0,1))

Could you take a look at this and see where I went wrong? Your help will be greatly appreciated,
Yaoping Wang
PhD student in the Ohio State University

[helske]

There are multiply issues here.

1. In order to get same results as with arima function, you need to set H to zero, i.e. you do not want additional noise.
2. Your initial values to Q throws an error as you are starting from value exp(100)=2.688117e+43 which is way too large for any practical application (it would cause severe numerical problems, so model checking throws error without even attempting to run Kalman filter).
3. You need to restrict the AR coefficients to stationary region. See for example function artransform (note that I just found a memory leak there, I will upload new version to CRAN today).
4. You need to also update P1, the initial covariance matrix of the states.

Here's an example with Nile data, I omit the predifferencing and fit ARIMA(2,1,1) model:

data(Nile)
model <- SSModel(Nile ~ SSMarima(ar=c(0.1,0.2),ma=c(0.1),d=1), H=0)

#due to stationarity checks it's easier to use own objective function and optim directly
likfn <- function(pars, model, estimate=TRUE){
  # use artransform to transform parameters so the model is stationary and invertible
  tmp <- try(SSMarima(artransform(pars[1:2]),
                  artransform(pars[3]),d=1,Q = exp(pars[4])),silent=TRUE)    
  if(!inherits(tmp,"try-error")){
    model["T","arima"] <- tmp$T 
    model["R","arima"] <- tmp$R    
    model["P1","arima"] <- tmp$P1
    model["Q","arima"] <- tmp$Q
    if(estimate){
      -logLik(model)
    } else model
  } else {
    if(estimate){
      1e100
    } else model
  } 
}

inits <- c(0.1,0.5,0.5,log(15000))
fit <- optim(inits, likfn, model=model, method='BFGS')
model <- likfn(fit$par,model,FALSE)

model$T
model$R
model$Q
#very similar results as with arima function:
arima(x=Nile, order = c(2,1,1))

Note that in some cases we get wildly different results due to different starting values, numerical instabilities near the stationary region, and because sigma^2 is estimated different way in arima (if I remember correctly it is computed from standardized residuals after running Kalman filter with unit variance).

[ypwong22]

Thank you and I updated my KFAS package. The example code worked for me, but there is a new question. When I change the initial parameters in the above example, e.g. inits <- c(0.1,0.5,0.5,log(15000)) changed to inits <- c(0.1,0.5,0.5,log(3000)), the estimated model parameters also changed greatly. Is there a reason why this happens? Is there a rule on prescribing the initial conditions?

[helske]

The initial value log(3000) is so far off that during the optimization routine the variance Q starts to go towards zero, which causes the model to become degenerate (as the H=0 as well). There is a check in logLik.SSModel against cases where all values of H and Q are smaller than machine epsilon, but it doesn't always help, like here where adding stricter check in likfn (like if(exp(pars[4])<1e-5) return(1e100)) still results a case where parameter estimates go towards the boundary.

In general state space models there isn't really rules for initial conditions, and the likelihood contains multiple maximums, so basically you try to guess something and test different set-ups and see what gives you best likelihood (and that the values are not insensible like in above example). For simple models you can try to infer the initial values from the data, for example here you have just one variance component, so good initial value for that could be log(var(Nile)). For ARMA parameters, zero is often a good first try, or some other small values around zero (using artransform for checking the values after transformation).

[helske]

I guess this can be closed now.