Extrapolating covariate for hierarchical model

[andrewbcooper]

Hello again! I am trying to model a variable (y) as a function of a trend and a covariate (x). However, I do not have values of the covariate for all future predictions of y. I've come up with brute-force way of predicting future Xs, and then passing those simulations to multivariate model to get future predictions of Y. However, this is very slow and not very elegant. I'm thinking there must be a more straight-forward way to do this. (Note, I'm doing it all in log-space in order to capture geometric growth.) Do you have any suggestions? I've attached sample code below. Thanks! Cheers, Andy

library(reshape2)
library(KFAS)

MyData <- data.frame(Year = 2011:2020,
MyY = c(3437,4422,4665,4983,5459,5925,NA,NA,NA,NA),
MyX = c(2242,3033,3158,3348,3602,4102,4135,NA,NA,NA))

Build the model for X

MyModelX <- SSModel((log(MyData$MyX)) ~
SSMtrend(degree=2,Q = list(matrix(NA),matrix(NA))), H=NA)

Fit the model for X

MyFitX <- fitSSM(MyModelX, inits = c(rep(0, 3)),method="L-BFGS-B")

Simulate the singal for X

MySimsSignalX <- simulateSSM(MyFitX$model, type = "signals", nsim = 100,
antithetics = FALSE)

Simulate the epsilons for X

MySimsEpsilonX <- simulateSSM(MyFitX$model, type = "epsilon", nsim = 100,
antithetics = FALSE, conditional = FALSE)

Calculate the future values for X

MySimsX <- data.frame((MySimsSignalX[ , 1, ] + MySimsEpsilonX[ , 1, ]))

Replace simulated past values with actual observed values of X

MySimsX[1:7,] <- log(MyData$MyX[1:7])

Reformat data for multivariate KFAS

datay <- matrix(log(MyData$MyY),nrow=length(MyData$MyY),ncol=100)
meltedX <- melt(MySimsX)
datax <- split(meltedX,meltedX$variable)

Build model for trend in Y along with regression on log(X)

MyModelY <- SSModel(datay ~-1+SSMregression(rep(list(~value),100),Q=matrix(NA,1,1),type="common",data=datax) +
SSMtrend(degree=2,type="common",Q=c(NA,NA)),
H=diag(NA,100))

Fit the model

MyFitY <- fitSSM(MyModelY, inits = c(rep(0, 105)),method="L-BFGS-B")

Simulate the signal for Y

MySimsSignalY <- simulateSSM(MyFitY$model, type = "signals", nsim = 50,
antithetics = FALSE)

Simulate the epsilons for Y

MySimsEpsilonY <- simulateSSM(MyFitY$model, type = "epsilon", nsim = 50,
antithetics = FALSE, conditional = FALSE)

Calculate the forecasted Ys

MySimsY <- data.frame(exp(MySimsSignalY + MySimsEpsilonY))

[helske]

First possibility that comes to my mind would be the use of so called "seemingly unrelated time series model", where you use multivariate model for both Y and X, i.e. something like this:

model <- SSModel(log(cbind(MyY, MyX)) ~ 
  SSMtrend(2, Q = list(matrix(NA, 2, 2), matrix(NA, 2, 2))), 
  H = matrix(NA, 2, 2), data = MyData)

fit <- fitSSM(model, inits = rep(0, 9), method = "BFGS")
sim <- simulateSSM(fit$model, "obs", nsim = 100)
ts.plot(exp(sim[,1,]))

Essentially you are modelling both the observations Y and the covariate values X, instead of treating X fixed, and here the relationship between Y and X is incorporated into the covariance terms in H and Q.

[andrewbcooper]

Brilliant! It never occurred to me to incorporate the relationship between X and Y in the covariance terms. And yes, treating those future Xs as "known" in the model for Y was bothering me, but I couldn't figure out how to avoid it. This is perfect! And it takes 1/60th of the processing time! Thanks!