Question building an Ar(2) model

[shankykm]

Hello there, I am currently using KFAS package to build an AR(2) model on a univariate dummy data. I am building this model using two different methods, SSMarima and SSMcustom function, to test my understanding. However, the AR coefficients value from the two methods are not matching. From the output, you can see that the param values from SSMcustom method is close to the actual simulated values but not the SSMarima values. I am still learning SSM and KFAS package. So, my knowledge is rusty. Can someone please help me out? Below is the R code I am using. Thank you and appreciate any guidance/feedback.

################### Kalman filter Univariate AR(2) model ################

# simulate random data
n <- 100
Y <- arima.sim(n = n, model = list(ar = c(0.6, -0.4)))

# Ht = 1  # Measurement error variance
Zt = matrix(c(1,0), ncol = 2)  # Measurement matrix
Tt = matrix(c(0, 1, 0, 0), ncol = 2)  # AR(2) state equation matrix
# Rt <- diag(2)*0.01
Qt <- 1 #diag(2)*0  # State innovation covariance
a1 = rep(0, 2)  # Initial state mean
P1 = diag(2)  # Initial state covariance

# Method 1: Define an AR(2) univariate state space model using SSMcustom function
model_ar2_1 <- SSModel(Y ~ -1 + 
                      SSMcustom(Z = Zt, T = Tt, Q=Qt))

# Method 2: Define an AR(2) univariate state space model using SSMarima function
model_ar2_2 <- SSModel(Y ~ -1 + 
                      SSMarima(ar=c(0, 0), d=0,  Q=1))

# Update function for method 1
updatefn_1 <- function(pars, model){
  print(pars)
  # model$H[1] <- pars[1] # measurement noise
  model$T[, , 1][1,1] <- pars[1] # AR1 coefficient
  model$T[, , 1][1,2] <- pars[2] # AR2 coefficient
  model$Q[, , 1][1] <- pars[3] # process noise
  model
}

# Update function for method 2
updatefn_2 <- function(pars, model){
  print(pars)
  tmp <- SSMarima(ar = pars[1:2], d = 0, Q = pars[3])
  model["R", states = "arima"] <- tmp$R
  model["Q", states = "arima"] <- tmp$Q
  model["P1", states = "arima"] <- tmp$P1
  model
}

# check that variances are non-negative
check_model <- function(model) {
  (model["H"] > 0 && model["Q"] > 0)
}

# Method 1 fitSSM
fit_ar2_1 <- fitSSM(model_ar2_1, c(0.01, 0.01, 1), updatefn_1,
                   lower=c(-0.99, -0.99, 1), upper=c(0.99, 0.99, 100), 
                   method = "L-BFGS-B", checkfn = check_model)

# Method 2 fitSSM
fit_ar2_2 <- fitSSM(model_ar2_2, c(0.01, 0.01, 1), updatefn_2,
                    lower=c(-0.99, -0.99, 1), upper=c(0.99, 0.99, 100), 
                    method = "L-BFGS-B", checkfn = check_model)

#Method 1 Parm output
fit_ar2_1$optim.out$par

#Method 2 Parm output
fit_ar2_2$optim.out$par

#Method 1 final state estimates
(out_arima_1 <- KFS(fit_ar2_1$model))

#Method 2 final state estimates
(out_arima_2 <- KFS(fit_ar2_2$model))

[helske]

You are not updating T in your second update function.

[shankykm]

Thanks, Dr. Helske. That did work. I have one follow up question. How should my second update function look like if i were to specify a non constant variance, Q, for my time series? For example, in my above simulated data (Y), if I want to specify that 1:40, 41:70 and 71:100 time points should have different variance, Q, then how could i incorporate this in my second update function for SSMarima function? Thanks, again.

[helske]

The SSMarima function assumes time-invariant Q, so you should use SSMcustom and then in the update function do something like

 model$Q[, , 1:40] <- pars[3]
 model$Q[, , 41:70] <- pars[4]
 model$Q[, , 71:100] <- pars[5]

where the original model is built with Q being 1x1x100. You can also use SSMarima but then you need to modify not only Q but also tv attribute of the SSModel object, i.e. attr(model, "tv")[5] <- 1L (see ?SSModel), e.g., something like

model$Q <- array(x, c(1,1,100))
attr(model, "tv")[5] <- 1L

[shankykm]

Dr. Helske - Thanks for the response. That did work. As I was building my knowledge using the KFAS package, I stumbled upon another error and this time with bivariate AR(1) model. Below is the code I am using. What I want to do is estimate the AR coefficients separately for each series regressing on a specific covariate. I also wanted to estimate the time invariant covariance matrix Q for the two series. When i run this model I am getting "System matrices (excluding Z) contain NA or infinite values, covariance matrices contain values larger than 1e+07". I believe this is happening because my covariance update step in the update function is not occuring correctly and I don't know how to fix this. Can you please help me out so my model can estimate the two AR components, regressors and covariance matrix correctly? Thanks again for all your help.

# simulate random data
n <- 100
Y1 <- arima.sim(model = list(ar = .99)) # series 1
Y2 <- arima.sim(model = list(ar = .5)) # series 2

X1 <- rnorm(n, mean = 0, sd = 1) # covariate for series 1
X2 <- rnorm(n, mean = 0.5, sd = 1) # covariate for series 2

Qt <- matrix(rep(1, 4), ncol=2) #variance covariance matrix to be estimated for bivariate model. 

# Define an AR(2) bivariate state space model using SSMarima function with covariates
bivar_model_ar2 <- SSModel(cbind(Y1, Y2) ~ -1  +
                           SSMregression(~X1, index=c(1)) +
                           SSMregression(~X2, index=c(2)) +
                           SSMarima(ar=rep(0, 2), d=0, Q=Qt)
)

# Update function for model
updatefn <- function(pars, model){
  print(pars)
  tmp <- try(SSMarima(ar = artransform(pars[1:2]), d = 0,
             Q=exp(pars[3:6])), silent = TRUE)

  # stationary check, see note in artransform docs
  if(inherits(tmp, "try-error")) {
    model$Q[] <- NA # set something to NA just in case original model is ok
    return(model) # this goes to checkfn and causes rejection due to NA values
  }

  model["R", states = "arima"] <- tmp$R
  model["Q", states = "arima"] <- tmp$Q
  model["P1", states = "arima"] <- tmp$P1
  model["T", states = "arima"] <- tmp$T
  model
}

# fit the model to the dataset
fit_ar2 <- fitSSM(bivar_model_ar2, init=c(rep(0.01,2), rep(1, 4)),
                    method = "BFGS", updatefn=updatefn)

#Parm output
fit_ar2$optim.out$par

#final state estimates
(out_arima_2 <- KFS(fit_ar2$model))

[shankykm]

Dr. Helske - Just wanted to follow up and update on my previous question around estimating the parameters (AR components, regressors and covariance matrix) for a bivariate model. After reading from the R help and your documentation, I updated my below code. My questions are,

1. Is this the right way to estimate the parameters especially the AR(1) components and covariance matrix for a bivariate model?
2. if i were to predict future observations, do I need to provide the model for the future observations via newdata argument and supply all the necessary estimates to SSModel object? I did that in my below code and its throwing another error saying "Misspecified Z, argument Z must be a (p x m) matrix".

Thanks again for all your help.

# simulate random data
n <- 100
Y1 <- arima.sim(n, model = list(ar = .99)) # series 1
Y2 <- arima.sim(n, model = list(ar = -.5)) # series 2

plot(Y1)
plot(Y2)

# covariates
X1 <- rnorm(n, mean = 0, sd = 100) # covariate for series 1
X2 <- rnorm(n, mean = 0.5, sd = 100) # covariate for series 2
Y = cbind(Y1, Y2)

Qt <- matrix(rep(1, 4), ncol=2) #variance covariance matrix to be estimated for bivariate model. 

# Define an AR(2) bivariate state space model using SSMarima function with covariates
bivar_model_ar2 <- SSModel(Y ~ -1  +
                           SSMregression(~X1, index=c(1), Q=NA) +
                           SSMregression(~X2, index=c(2), Q=NA) +
                           SSMarima(ar=rep(0, 1), d=0, 
                                    Q=matrix(rep(NA, 4), ncol=2)
                                    )
                           )

# Update function for model
updatefn <- function(pars, model){
  T <- artransform(pars[1:2])
  model["T", states = "arima"] <- T
  Q <- diag(pars[3:4])
  Q[upper.tri(Q)] <- pars[5]
  model["Q", etas = "arima"] <- model["P1", states = "arima"] <- crossprod(Q)
  Q <- diag(pars[6:7])
  model["Q", etas = "regression"] <- crossprod(Q)
  model["R", ] <- model$R
  model
}

# check that variances are non-negative
check_model <- function(model) {
  (diag(model["Q"][,,1]) > 0 && model["H"] > 0)
}

# fit the model to the dataset
fitinit <- fitSSM(bivar_model_ar2, init=c(rep(0.01,2), rep(5, 5)),
                  lower=c(rep(-0.99, 2), rep(1, 5)), 
                  upper=c(rep(0.99, 2), rep(10, 5)),
                  method = "L-BFGS-B", 
                  updatefn=updatefn, checkfn = check_model)

fitinit$optim.out$par

#refit model using the estimated params as initial state values
fit <- fitSSM(bivar_model_ar2, init=fitinit$optim.out$par,
              nsim = 250,
              method = "BFGS", 
              updatefn=updatefn, checkfn = check_model)

#Parm output
fit$optim.out$par

# final state estimates
(out_arima_2 <- KFS(fit$model, nsim = 1000))

# var covar matrix estimates
(varcor <- fit$model["Q", etas = "arima"])
(varcor <- fit$model["Q", etas = "regression"])

# new prediction
pred <- predict(fit$model, 
                newdata = SSModel(ts(matrix(NA, 6, 2)) ~ -1 +
                                  SSMregression(~X1[1:6], index=c(1)) +
                                  SSMregression(~X2[1:6], index=c(2)) +
                                  SSMcustom(Z = fit$model$Z, 
                                            T = fit$model$T, 
                                            R = fit$model$R,
                                            Q = fit$model$Q)
                                  ),
                interval = "confidence", nsim = 10000)

# ggplot of actual vs predict
cbind(data.frame(predict(fit$model, interval="prediction", level=0.95)), Y1=Y1, Y2=Y2, time=1:n) %>% 
  ggplot(., aes(x=time, y=Y1))+
   geom_point(aes(colour="Y1"), size=1) +
   geom_line(aes(x=time, y=Y1.fit, colour="Predict"), size=0.8) +
  geom_ribbon(aes(ymin = Y1.lwr, ymax = Y1.upr), alpha = 0.4)

cbind(data.frame(predict(fit$model, interval="prediction", level=0.95)), Y1=Y1, Y2=Y2, time=1:n) %>% 
  ggplot(., aes(x=time, y=Y2))+
  geom_point(aes(colour="Y1"), size=1) +
  geom_line(aes(x=time, y=Y2.fit, colour="Predict"), size=0.8) +
  geom_ribbon(aes(ymin = Y2.lwr, ymax = Y2.upr), alpha = 0.4)

[helske]

In updatefn you are incorrectly updating the T matrix: That contains four elements but you are only assigning two. model["T", states = "arima"] <- diag(T) is one way to fix the issue. Also the update on R is not necessary here, and the second run of optimisation is also not necessary as argument nsim is ignored in case of linear-Gaussian model like yours (it also has no effect in predict).

In predict the SSMcustom is not correct as fit$model$Z etc contains not only the SSMarima part but the whole model so the old model and the new model are not compatible; you should build the model for newdata same way as you did for the orignal model. Or, another option build a new model with missing observations in the end corresponding to those time points you want to predict and use that in fitting and subsequent prediction without newdata argument.

PS. The github issues are mainly for reporting bugs or feature requets, this type of general questions regarding the usage of the package are more suited for example to stackoverflow (you can email me with link to SO question if it looks like you are not getting answers there).