Esempio stima parametri ignoti con KFAS

Articoli di approfondimento

Problema

Per questra applicazione creiamo noi una serie storica:

\[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.
# 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))

model1$optim.out$convergence
## [1] 0
Q_vera<-matrix(c(sigma2_nu,0,0,0,sigma2_zeta,0,0,0,sigma2_omega),ncol=3)
cat("\n \n  Matrice Q \n");Q_vera;cat("\n \n Matrice Q prevista \n");model1$model$Q
## 
##  
##   Matrice Q
##           [,1]      [,2]         [,3]
## [1,] 0.1353353 0.0000000 0.000000e+00
## [2,] 0.0000000 0.2231302 0.000000e+00
## [3,] 0.0000000 0.0000000 4.539993e-05
## 
##  
##  Matrice Q prevista
## , , 1
## 
##           [,1]      [,2]         [,3]
## [1,] 0.1510326 0.0000000 0.000000e+00
## [2,] 0.0000000 0.2093234 0.000000e+00
## [3,] 0.0000000 0.0000000 1.198907e-07
cat("\n \n  Matrice H \n");sigma2_epsilon;cat("\n \n Matrice H prevista \n");model1$model$H
## 
##  
##   Matrice H
## [1] 1.491825
## 
##  
##  Matrice H prevista
## , , 1
## 
##          [,1]
## [1,] 1.597679
  1. Generare i residui ausiliari per costruire i \(t\)-test per identificare outlier additivi (usare la funzione rstandard(oggettoKFS, "pearson"))
fit1 <- KFS(model1$model)
residui <- rstandard(fit1, "pearson")
plot(residui)
abline(h = c(sort(residui)[as.integer(0.975*length(residui))], sort(residui)[as.integer(0.025*length(residui))]), lty = 2, col = "red", lwd = 2)

  1. Stimare un modello UCM (model2) per il quale alle componenti di cui sopra si aggiunga una variabile dummy per ogni residuo ausiliare maggiore di 2.8 in valore assoluto.
outlier <- which(abs(residui)>2.8)
dum <- matrix(0, nrow=length(y), ncol=length(outlier))
for (j in 1:length(outlier)) dum[outlier[j], j] <- 1

mod2 <- SSModel(y~0+dum+SSMtrend(2, list(NA, NA))+SSMseasonal(4, NA), H = NA)
model2 <- fitSSM(mod2, c(sigma2_nu, sigma2_zeta, sigma2_omega, sigma2_epsilon))

model2$optim.out$convergence
## [1] 0
cat("\n \n  Matrice Q \n");Q_vera;cat("\n \n Matrice Q prevista \n");model2$model$Q
## 
##  
##   Matrice Q
##           [,1]      [,2]         [,3]
## [1,] 0.1353353 0.0000000 0.000000e+00
## [2,] 0.0000000 0.2231302 0.000000e+00
## [3,] 0.0000000 0.0000000 4.539993e-05
## 
##  
##  Matrice Q prevista
## , , 1
## 
##           [,1]      [,2]       [,3]
## [1,] 0.1386115 0.0000000 0.00000000
## [2,] 0.0000000 0.2086863 0.00000000
## [3,] 0.0000000 0.0000000 0.00055488
cat("\n \n  Matrice H \n");sigma2_epsilon;cat("\n \n Matrice H prevista \n");model2$model$H
## 
##  
##   Matrice H
## [1] 1.491825
## 
##  
##  Matrice H prevista
## , , 1
## 
##          [,1]
## [1,] 1.464188
  1. Stimare un modello UCM (model3) identico al model1, ma dove la variabile dipendente è sostituita con una variabile identica, ma con NA che sostituiscono i valori in corrsipondenza dei quali i residui ausiliari sono maggiori di 2.8 in valore assoluto.
y_new <- y
y_new[outlier] <- NA
mod3 <- SSModel(y_new~0+SSMtrend(2, list(NA, NA))+SSMseasonal(4, NA), H = NA)
model3 <- fitSSM(mod3, c(sigma2_nu, sigma2_zeta, sigma2_omega, sigma2_epsilon))


model3$optim.out$convergence
## [1] 0
cat("\n \n  Matrice Q \n");Q_vera;cat("\n \n Matrice Q prevista \n");model3$model$Q
## 
##  
##   Matrice Q
##           [,1]      [,2]         [,3]
## [1,] 0.1353353 0.0000000 0.000000e+00
## [2,] 0.0000000 0.2231302 0.000000e+00
## [3,] 0.0000000 0.0000000 4.539993e-05
## 
##  
##  Matrice Q prevista
## , , 1
## 
##           [,1]      [,2]       [,3]
## [1,] 0.1386115 0.0000000 0.00000000
## [2,] 0.0000000 0.2086863 0.00000000
## [3,] 0.0000000 0.0000000 0.00055488
cat("\n \n  Matrice H \n");sigma2_epsilon;cat("\n \n Matrice H prevista \n");model3$model$H
## 
##  
##   Matrice H
## [1] 1.491825
## 
##  
##  Matrice H prevista
## , , 1
## 
##          [,1]
## [1,] 1.464188