Application: Estimating Unknown Parameters with KFAS

In-depth Articles

Problem

For this application we create a time series:

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

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 the residuals ausiliari per costruire the \(t\)-test per identificare outlier additivi (usare the function rstandard(oggettoKFS, "pearson"))
fit1 <- KFS(model1$model)
residuals <- rstandard(fit1, "pearson")
plot(residuals)
abline(h = c(sort(residuals)[as.integer(0.975*length(residuals))], sort(residuals)[as.integer(0.025*length(residuals))]), lty = 2, col = "red", lwd = 2)

  1. Stimare a model UCM (model2) per the quale to the components di cui sopra it aggiunga a variable dummy for every residuo ausiliare maggiore di 2.8 in value assoluto.
outlier <- which(abs(residuals)>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 a model UCM (model3) identico al model1, ma where the variable dependent is sostituita with a variable identica, ma with NA che sostituiscono the values in corrsipondenza of the quali the residuals ausiliari are maggiori di 2.8 in value 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