Forecasting, filtering e smoothing

Una volta stimati i parametri ignoti del modello in forma state space, è naturale volere usare il modello per prevedere le serie storiche o fare inferenza sulle loro componenti.

Per questra applicazione creiamo noi una serie storica con UCM cosi da poter controllare i risultati:

\[Y_t= \mu_t+\gamma_t+ \epsilon_t \\.\\ Dove: \\ \mu_t = \mu_{t-1}+\beta_t+ \nu_t \space \space con \space \space \nu_t \sim WN(0, \sigma^2_{\nu})\\ \beta_t = \beta_{t-1}+ \zeta_t \space \space con \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 con \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)
# creiamo una serie storica mensile per 20 anni

20*12

## [1] 240

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

#costruzione sinusoide: cos(freq*tempo) & sen(freq*tempo)
stagionesinus<-cbind(cos(outer(1:240,frequenze_stagionali)),sin(outer(1:240,frequenze_stagionali)))
#l'ultima colonna va tolta perchè è sempre 0



beta_t<-mu_t<-0
for(i in 2:240){
  
  set.seed(as.integer(15121/i*101)) # per riproducibilita
  beta_t[i]<-beta_t[i-1]+ rnorm(1,mean=0, sd=sqrt(sigma2_zeta))
  
  set.seed(as.integer(15121/i*13)) # per riproducibilita
  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)

Suggerimento: come valori iniziali per il log delle varianze usate \[ \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. Stimare un modello UCM (model1) composto da local linear trend, stagionalità a dummy stocastiche, errore di osservazione. (esercizio fatto nei capitoli precedenti: Link dell'articolo)

# sono tutti in logaritmo
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)
residui <- rstandard(fit1, "pearson")
outlier <- which(abs(residui)>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))

Ipotiziamo una serie storica: \(Y_T= \{y_1,y_2,...,y_t\}\) In questo ambito l’inferenza sulle variabili di stato prende nomi diversi in base all’informazione che si utilizza:

si parla di filtro (filtering; previsione in tempo reale) se si calcola:

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

si parla di lisciatura (smoothing) se si calcola:

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

si parla di previsione (forecasting) se si calcola:

\[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)) # Grafico della serie storica 
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)) # Grafico della serie storica 
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)) # Grafico della serie storica 
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)) # Grafico della serie storica 
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)) # Grafico della serie storica 
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)) # Grafico della serie storica 
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="Previsioni")
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