Forecasting, Filtering and Smoothing

Once the unknown parameters of the state space model have been estimated, it is natural to want to use the model to forecast the time series or make inference on their components.

For this application we create a time series with a UCM so that we can verify the results:

\[Y_t= \mu_t+\gamma_t+ \epsilon_t \\.\\ Where: \\ \mu_t = \mu_{t-1}+\beta_t+ \nu_t \space \space with \space \space \nu_t \sim WN(0, \sigma^2_{\nu})\\ \beta_t = \beta_{t-1}+ \zeta_t \space \space with \space \space \zeta_t \sim WN(0, \sigma^2_{\zeta})\\ \gamma_t = \sum_{j=1}^{2}{20 \cos{\biggl(\frac{2 \pi}{s} \space j \space t\biggr)}+10\sin{\biggl(\frac{2 \pi}{s} \space j \space t\biggr)}}+ \omega_t \space \space with \space \space \omega_t \sim WN(0, \sigma^2_{\omega}) \\ \epsilon_t = WN(0, \sigma^2_{\epsilon}) \]

log_sigma2_nu=-2
log_sigma2_zeta=-1.5
log_sigma2_omega=-10
log_sigma2_epsilon=0.4

sigma2_nu=exp(-2)
sigma2_zeta=exp(-1.5)
sigma2_omega=exp(-10)
sigma2_epsilon=exp(0.4)
# we create a monthly time series for 20 years

20*12

## [1] 240

mu_t<-1:240
frequencies_seasonal<-2*pi*(1:2)/4

# sinusoid construction: cos(freq*time) & sin(freq*time)
stagionesinus<-cbind(cos(outer(1:240,frequencies_seasonal)),sin(outer(1:240,frequencies_seasonal)))
# the last column must be removed because it is always 0



beta_t<-mu_t<-0
for(i in 2:240){
  
  set.seed(as.integer(15121/i*101)) # for reproducibility
  beta_t[i]<-beta_t[i-1]+ rnorm(1,mean=0, sd=sqrt(sigma2_zeta))
  
  set.seed(as.integer(15121/i*13)) # for reproducibility
  mu_t[i]<-mu_t[i-1]+beta_t[i]+ rnorm(1,mean=0, sd=sqrt(sigma2_nu))
  
  
}
set.seed(15121)
gamma_t<-stagionesinus %*% c(rep(20,2),rep(10,2))+rnorm(240,mean=0, sd=sqrt(sigma2_omega))

set.seed(15121*11)
y<-mu_t+gamma_t+rnorm(240,sd=sqrt(sigma2_epsilon))
ts.plot(y)

Hint: as initial values for the log of the variances use \[ \text{(level)} \log\sigma^2_\eta = 0, \;\; \text{(slope)} \log\sigma^2_\zeta = -1.5, \;\; \text{(seas)}\log\sigma^2_\omega = 5, \;\; \text{(irregular)} \log\sigma^2_\epsilon = 0.4 \]

1. Estimate a UCM model (model1) composed of a local linear trend, stochastic dummy seasonality, and observation error. (exercise from the previous chapters: Article link)

# are all in logarithm
log_sigma2_nu=0
log_sigma2_zeta=-1.5
log_sigma2_omega=0
log_sigma2_epsilon=0.4
library(KFAS)
mod1 <- SSModel(y~0+SSMtrend(2, list(NA, NA))+SSMseasonal(4, NA), H = NA)
model1 <- fitSSM(mod1, c(log_sigma2_nu, log_sigma2_zeta, log_sigma2_omega, log_sigma2_epsilon))
fit1 <- KFS(model1$model)
residuals <- rstandard(fit1, "pearson")
outlier <- which(abs(residuals)>2.8)
for(i in 1:length(outlier)) y[outlier[i]] <- mean(c(y[outlier[i]-1],y[outlier[i]+1]))
mod1 <- SSModel(y~0+SSMtrend(2, list(NA, NA))+SSMseasonal(4, NA), H = NA)
model1 <- fitSSM(mod1, c(sigma2_nu, sigma2_zeta, sigma2_omega, sigma2_epsilon))

Let us consider a time series: \(Y_T= \{y_1,y_2,...,y_t\}\) In this ambito l’inferenza on the variables di stato prende nomi diversi in base all’informazione che it utilizza:

We speak of filtering (real-time prediction) when we compute:

\[a_{t|t} = P[a_t|Y_t] \]

We speak of smoothing when we compute:

\[a_{t|s} = P[a_t|Y_s] \space \space t.c. \space \space t<s\]

We speak of forecasting when we compute:

\[a_{t|s} = P[a_t|Y_s] \space \space t.c. \space \space t>s\]

filt1 <- KFS(model1$model, filtering = c("state", "signal"))
fitted<-data.frame(filt1$alphahat)

# level
ts.plot(mu_t, ylab=expression(mu)) # Plot of the time series 
lines(filt1$alphahat[,"level"], col = "blue")

print(mean(sum((mu_t-fitted$level)^2)))

## [1] 240.931

# slope
ts.plot(beta_t, ylab=expression(beta)) # Plot of the time series 
lines(fitted$slope, col = "blue")

print(mean(sum((beta_t-fitted$slope)^2)))

## [1] 61.51514

# stag
ts.plot(gamma_t, ylim=c(-50,50), ylab=expression(gamma)) # Plot of the time series 
lines(rowSums(fitted[,3:5]), col = "coral") 

print(mean(sum((gamma_t-rowSums(fitted[,3:5]))^2)))

## [1] 120162

smo1 <- KFS(model1$model, smoothing = c("state","signal")) 
fitted<-data.frame(smo1$alphahat)



# level
ts.plot(mu_t, ylab=expression(mu)) # Plot of the time series 
lines(smo1$alphahat[,"level"], col = "blue")

print(mean(sum((mu_t-fitted$level)^2)))

## [1] 240.931

# slope
ts.plot(beta_t, ylab=expression(beta)) # Plot of the time series 
lines(fitted$slope, col = "blue")

print(mean(sum((beta_t-fitted$slope)^2)))

## [1] 61.51514

# stag
ts.plot(gamma_t, ylim=c(-50,50), ylab=expression(gamma)) # Plot of the time series 
lines(rowSums(fitted[,3:5]), col = "coral") 

print(mean(sum((gamma_t-rowSums(fitted[,3:5]))^2)))

## [1] 120162

train<-y[1:180]
test<-ts(y[181:240],start = 181)


mod2 <- SSModel(train~0+SSMtrend(2, list(NA, NA))+SSMseasonal(4, NA), H = NA)
model2 <- fitSSM(mod2, c(log_sigma2_nu, log_sigma2_zeta, log_sigma2_omega, log_sigma2_epsilon))




prev <- predict(model2$model, n.ahead = length(test), interval = "prediction")

ts.plot(prev, col=c(1,2,2), main="Forecasts")
lines(test,col=4)

ts.plot(test-prev[,1], col=c(1,2,2), main="Errori")
abline(h=0,col="green")

MSE_tot<-mean(sum((test-prev[,1])^2));message(MSE_tot)

## 8765.85476440955