Modelli UCM in forma state space in KFAS

Articoli di approfondimento

Modelli state space in KFAS

Definizione generale

La forma state space utilizzata da KFAS è leggermente diversa da quella che abbiamo usato fino ad ora. In particolare non è possibile utilizzare i vettori ct e dt e la matrice di covarianza dei disturbi del vettore di stato è specificata utilizzando il prodotto di due matrici diverse. KFAS usa le seguente forma state space:

  • Equazioni di misurazione (o di osservazione): mettono in relazione lineare una serie storica (o un vettore di serie storiche) con \(m\) variabili di stato che contengono le componenti non osservabili nel caso degli UCM e altre componenti che aiutano a costruirle.\[ \underline{Y}_t = \underline{Z}_t \space \underline{\alpha}_t + \underline{\epsilon}_t \space, \space \space con \space \space \underline{\epsilon}_t \sim WN(\underline{0}, H_t) \] con \(\alpha_t\) vettore di stato che non deve avere necessariamente la dimensione della serie storica.
  • Equazioni di stato (o di transizione): informano circa la legge del moto, cioè come \(\alpha_t\) evolve nel tempo (si tratta di un VAR(1) con coefficienti che possono evolvere nel tempo).\[ \underline{\alpha}_{t+1} = \underline{T}_t \space \underline{\alpha}_t + \underline{\nu}_t \space, \space \space con \space \space \underline{\nu}_t = \underline{R}_t \underline{\eta}_t \sim WN(\underline{0}, Q_t) \]Tale forma è chiamata forma futura perchè a sinistra c’è solo il tempo \(t+1\) e a destra il tempo \(t\). Solitamente si scrive \(\ni_{t+1}\) ma trattandosi di un WN non cambia nulla. Il compito della matrice \(R\) è quello di attribuire i due disturbi alle equazioni di transizioni giuste: ad esempio associare il primo disturbo alla seconda equazione e il secondo disturbo alla terza equazione.

La forma state-space è completata dalla specificazione dei primi due momenti di \(\alpha_1\): \[ \underline{a}_{1|0} = E(\underline\alpha_1) \\ P_{1|0}=E[(\underline\alpha_1-\underline\alpha_{1|0})(\underline\alpha_1-\underline\alpha_{1|0})^T \space]\]

Nel caso di componenti stazionarie si inseriscono la media e la varianza di tali componenti, nel caso invece di componenti non stazionarie (come trend o la stagionalità) si possono inserire condizioni diffuse, cioè una media arbitraria e delle varianze infinite (c’è infinita incertezza rispetto al valore che la v.c. può assumere). La possibilità di porre condizioni diffuse garantisce che definire \(\alpha_{1|0}\) e \(P_{1|0}\) non sia un limite ma anche l’algoritmo usato non lo consente si può fissare una varianza così alta da esse indistinguibile da infinito da punto di vista pratico.

Condizioni tecniche:

  • Prima condizione: \(E(\underline\epsilon_t \space \underline\nu_t^T)=G_t\) e si suppone che \(G_t=0\) per semplicità;
  • Seconda condizione: \(E(\underline\epsilon_t \space \underline\nu_s^T)=0\) per ogni \(t \neq s\);
  • Terza condizione: \(E[(\underline\alpha_1-\underline\alpha_{1|0})\space \underline\epsilon_t^T]=0\);
  • Quarta condizione: \(E[(\underline\alpha_1-\underline\alpha_{1|0})\space \underline\nu_t^T]=0\).

Definizione UCM in KFAS

Il pacchetto R che utilizzeremo nella stima di modelli in forma state space è KFAS, che potete scaricare dal CRAN come tutti i principali pacchetti R. Il pacchetto contiene degli ottimi algoritmi per fare inferenza sulle variabili di stato (filtro di Kalman e smoothers). Il pacchetto KFAS definisce una classe di oggetti di nome SSModel, che di fatto è una lista contenente tutte le informazioni della forma state space. Il costruttore, cioè la funzione che genera oggetti di tipo SSModel, si chiama a sua volta SSModel() e, ai nostri fini, necessita di soli tre parametri formali: formula, data e H. Quindi, per generare una forma state space e assegnarla alla variabile modello scriveremo:

modello <- SSModel(formula, data, H)

Il secondo e il terzo parametro formale sono molto banali: data deve contenere un data.frame con i dati necessarialmodello,cioè la o le variabili dipendenti e gli eventuali regressori(se le variabili sono disponibili nel global environment il parametro formale data può essere omesso); H deve contenere la matrice di covarianza (o l’array di matrici di covarianza se queste evolvono nel tempo) dell’errore di osservazione. Al parametro formale formula deve essere passato un oggetto di tipo formula comune in R per gestire variabili esogene ed endogene. Per generare una forma state space partendo direttamente dalle matrici di sistema si deve utilizzare all’interno del parametro formula la funzione:

SSMcustom(Z, T, R, Q, a1, P1, P1inf)

In KFAS esistono diverse funzioni utilizzabili additivamente per creare le componenti più comuni dei modelli a componenti non osservabili. Per i dettagli delle funzioni si consiglia di riferirsi alla documentazione del pacchetto KFAS, qui ci limitiamo a mostrare degli esempi per l’utilizzo in modello univariati. SSMarima(ar, ma, d = 0, Q) produce la forma state space di un modello ARIMA(p,d,q) con i coefficienti della parte AR e MA specificati nei vettori ar e ma, rispettivamente.

Nella prossima parte vedremo nel dettaglio la specificazione delle varie componenti dei modelli UCM nello specifico.

RW + noise

\[T=1, \space \space Q=\sigma_\eta^2, \space \space H=\sigma_\epsilon^2, \space \space Z=1 \] Quindi il modello state space diventa:

  • Equazioni di misurazione (o di osservazione):

\[ Y_t = \alpha_t + \epsilon_t \space, \space \space con \space \space \epsilon_t \sim WN(0, \sigma_\epsilon^2) \] - Equazioni di stato (o di transizione):

\[ \alpha_{t+1} = \alpha_t + \nu_t \space, \space \space con \space \space \nu_t \sim WN(0, \sigma_\eta^2) \]

library(KFAS)
y<-AirPassengers
sigma2_eta=10
sigma2_epsilon=2

modelloRW<-SSModel(y~0+SSMtrend(degree = 1,Q=sigma2_eta), H= sigma2_epsilon)


modelloRW$T
## , , 1
## 
##       level
## level     1
modelloRW$Z
## , , 1
## 
##      level
## [1,]     1
modelloRW$H
## , , 1
## 
##      [,1]
## [1,]    2
modelloRW$Q
## , , 1
## 
##      [,1]
## [1,]   10

RWD + noise (local linear trend)

\[T=\bigg[\begin{array}{cc} 1 &1 \\0&1\end{array} \bigg], \space \space Q=\bigg[\begin{array}{cc} \sigma_\eta^2&1 \\0&1\end{array} \bigg], \space \space H=1, \space \space Z=\big[\begin{array}{cc} 0&1\end{array} \big] \] Quindi il modello state space diventa:

  • Equazioni di misurazione (o di osservazione):

\[ Y_t = \big[\begin{array}{cc} 0&1\end{array} \big] \bigg[\begin{array}{cc} \alpha_{1,t+1} \\ \alpha_{2,t+1} \end{array} \bigg] + \epsilon_t \space, \space \space con \space \space \epsilon_t \sim WN(0, 1) \] - Equazioni di stato (o di transizione):

\[ \bigg[\begin{array}{cc} \alpha_{1,t+1} \\ \alpha_{2,t+1} \end{array} \bigg] = \bigg[\begin{array}{cc} \alpha_{1,t} \\ \alpha_{2,t}\end{array} \bigg]+ \bigg[\begin{array}{cc} \nu_t \\0 \end{array}\bigg] \space, \space \space con \space \space \nu_t \sim WN(0, \sigma_\eta^2) \]

library(KFAS)
y<-AirPassengers
sigma2_eta=10
sigma2_epsilon=2

modelloRWD<-SSModel(y~0+SSMtrend(degree = 2,Q=list(sigma2_eta,0)), H= 1)


modelloRWD$T
## , , 1
## 
##       level slope
## level     1     1
## slope     0     1
modelloRWD$Z
## , , 1
## 
##      level slope
## [1,]     1     0
modelloRWD$H
## , , 1
## 
##      [,1]
## [1,]    1
modelloRWD$Q
## , , 1
## 
##      [,1] [,2]
## [1,]   10    0
## [2,]    0    0

RW + ciclo

\[T=\begin{bmatrix}1 &0&0 \\0& \cos \frac{2 \space \pi}{periodo} & \sin \frac{2 \space \pi}{periodo} \\0 & -\sin \frac{2 \space \pi}{periodo} & \cos \frac{2 \space \pi}{periodo}\end{bmatrix} , \space \space Q=\begin{bmatrix} \sigma_\eta^2&0&0\\0& \sigma_k^2&0\\0&0& \sigma_k^2 \end{bmatrix} , \space \space H=0, \space \space Z=\big[\begin{array}{cc} 1&1&0\end{array} \big] \] Quindi il modello state space diventa:

  • Equazioni di misurazione (o di osservazione):

\[ Y_t = \big[\begin{array}{cc} 1&1&0\end{array}\big] \begin{bmatrix} \alpha_{t+1} \\ cycle_{t+1} \\ cycle_{t+1}^* \end{bmatrix} \] - Equazioni di stato (o di transizione):

\[ \begin{bmatrix} \alpha_{t+1} \\ cycle_{t+1} \\ cycle_{t+1}^* \end{bmatrix} = \begin{bmatrix}1 &0&0 \\0& \cos \frac{2 \space \pi}{periodo} & \sin \frac{2 \space \pi}{periodo} \\0 & -\sin \frac{2 \space \pi}{periodo} & \cos \frac{2 \space \pi}{periodo}\end{bmatrix} \begin{bmatrix} \alpha_{t} \\ cycle_{t} \\ cycle_{t}^* \end{bmatrix} +\begin{bmatrix} \nu_t \\k_t\\k_t \end{bmatrix} \space, \\ \space \space con \space \space \nu_t \sim WN(0, \sigma_\eta^2) \\ \space \space \space con \space \space k_t \sim WN(0, \sigma_k^2) \]

library(KFAS)
y<-AirPassengers
sigma2_eta=10
sigma2_epsilon=2
sigma2_k=0.1

periodo=11
modelloRWcycle<-SSModel(y~0+SSMtrend(degree = 1,Q=sigma2_eta)+SSMcycle(period =periodo,Q= sigma2_k ), H= 0)

modelloRWcycle$T
## , , 1
## 
##        level      cycle    cycle*
## level      1  0.0000000 0.0000000
## cycle      0  0.8412535 0.5406408
## cycle*     0 -0.5406408 0.8412535
modelloRWcycle$Z
## , , 1
## 
##      level cycle cycle*
## [1,]     1     1      0
modelloRWcycle$H
## , , 1
## 
##      [,1]
## [1,]    0
modelloRWcycle$Q
## , , 1
## 
##      [,1] [,2] [,3]
## [1,]   10  0.0  0.0
## [2,]    0  0.1  0.0
## [3,]    0  0.0  0.1

RWD + ciclo

\[T=\begin{bmatrix} 1 &1&0 &0 \\ 0&1&0&0 \\0&0& \cos \frac{2 \space \pi}{periodo} & \sin \frac{2 \space \pi}{periodo} \\0&0 & -\sin \frac{2 \space \pi}{periodo} & \cos \frac{2 \space \pi}{periodo}\end{bmatrix} , \space \space Q=\begin{bmatrix} \sigma_\eta^2&0&0&0\\0&0&0&0\\0&0& \sigma_k^2&0\\0&0&0& \sigma_k^2 \end{bmatrix} , \space \space H=0, \space \space Z=\big[\begin{array}{cc} 1&0&1&0\end{array} \big] \]

N.B. La matrice \(T\) contiene al suo interno la matrice di rotazione:\(R(y)\) Quindi il modello state space diventa:

  • Equazioni di misurazione (o di osservazione):

\[ Y_t = \big[\begin{array}{cc} 1&0&1&0\end{array}\big] \begin{bmatrix} \alpha_{1,t+1} \\ \alpha_{2,t+1} \\ cycle_{t+1} \\ cycle_{t+1}^* \end{bmatrix} \] - Equazioni di stato (o di transizione):

\[ \begin{bmatrix} \alpha_{1,t+1} \\ \alpha_{2,t+1}\\ cycle_{t+1} \\ cycle_{t+1}^* \end{bmatrix} = \begin{bmatrix} 1 &1&0 &0 \\ 0&1&0&0 \\0&0& \cos \frac{2 \space \pi}{periodo} & \sin \frac{2 \space \pi}{periodo} \\0&0 & -\sin \frac{2 \space \pi}{periodo} & \cos \frac{2 \space \pi}{periodo}\end{bmatrix} \begin{bmatrix} \alpha_{1,t} \\ \alpha_{2,t}\\ cycle_{t} \\ cycle_{t}^* \end{bmatrix} +\begin{bmatrix} \nu_t\\0 \\k_t\\k_t \end{bmatrix} \space, \\ \space \space con \space \space \nu_t \sim WN(0, \sigma_\eta^2) \\ \space \space \space con \space \space k_t \sim WN(0, \sigma_k^2) \]

library(KFAS)
y<-AirPassengers
sigma2_eta=10
sigma2_epsilon=2
sigma2_k=0.1

periodo=11
modelloRWDcycle<-SSModel(y~0+SSMtrend(degree = 2,Q=list(sigma2_eta,0))+SSMcycle(period =periodo,Q= sigma2_k ), H= 0)

modelloRWDcycle$T
## , , 1
## 
##        level slope      cycle    cycle*
## level      1     1  0.0000000 0.0000000
## slope      0     1  0.0000000 0.0000000
## cycle      0     0  0.8412535 0.5406408
## cycle*     0     0 -0.5406408 0.8412535
modelloRWDcycle$Z
## , , 1
## 
##      level slope cycle cycle*
## [1,]     1     0     1      0
modelloRWDcycle$H
## , , 1
## 
##      [,1]
## [1,]    0
modelloRWDcycle$Q
## , , 1
## 
##      [,1] [,2] [,3] [,4]
## [1,]   10    0  0.0  0.0
## [2,]    0    0  0.0  0.0
## [3,]    0    0  0.1  0.0
## [4,]    0    0  0.0  0.1

RW + stagionalità a dummy

Con le dummy vengono delle matrici impegnative perche lunghe quanto il numero dell dummy meno uno piu gli altri elementi del modello, quindi faremo solo l’esempio con 4 dummy stagionali e un random wlak cosi da semplificare la comprensione della componente stagionale nei modelli state space.

\[T=\begin{bmatrix} 1 &0&0 &0 \\0&-1& -1 & -1 \\0&1 & 0 & 0\\ 0&0&1&0\end{bmatrix} , \space \space Q=\begin{bmatrix} \sigma_\eta^2&0\\0&\sigma_\gamma^2 \end{bmatrix} , \space \space H=0, \space \space Z=\big[\begin{array}{cc} 1&1&0&0\end{array} \big] \]

N.B. La matrice \(T\) contiene al suo interno la matrice di rotazione:\(R(y)\) Quindi il modello state space diventa:

  • Equazioni di misurazione (o di osservazione):

\[ Y_t = \big[\begin{array}{cc} 1&1&0&0\end{array}\big] \begin{bmatrix} \alpha_{t+1} \\ dummy1_{t+1}\\ dummy2_{t+1} \\ dummy3_{t+1} \end{bmatrix} \] - Equazioni di stato (o di transizione):

\[ \begin{bmatrix} \alpha_{t+1} \\ dummy1_{t+1}\\ dummy2_{t+1} \\ dummy3_{t+1} \end{bmatrix} = \begin{bmatrix} 1 &0&0 &0 \\0&-1& -1 & -1 \\0&1 & 0 & 0\\ 0&0&1&0\end{bmatrix} \begin{bmatrix} \alpha_{t} \\ dummy1_{t}\\ dummy2_{t} \\ dummy3_{t} \end{bmatrix} +\begin{bmatrix} \nu_t\\\gamma_t \\\gamma_t\\\gamma_t \end{bmatrix} \space, \\ \space \space con \space \space \nu_t \sim WN(0, \sigma_\eta^2) \\ \space \space \space con \space \space \gamma_t \sim WN(0, \sigma_\gamma^2) \]

library(KFAS)
y<-AirPassengers
sigma2_eta=10
sigma2_epsilon=2
sigma2_k=0.1
sigma2_gamma=1

periodo=4
modelloRWdummy<-SSModel(y~0+SSMtrend(degree = 1,Q=sigma2_eta)+SSMseasonal(period =periodo,Q= sigma2_gamma, sea.type = "dummy" ), H= 0)

modelloRWdummy$T
## , , 1
## 
##            level sea_dummy1 sea_dummy2 sea_dummy3
## level          1          0          0          0
## sea_dummy1     0         -1         -1         -1
## sea_dummy2     0          1          0          0
## sea_dummy3     0          0          1          0
modelloRWdummy$Z
## , , 1
## 
##      level sea_dummy1 sea_dummy2 sea_dummy3
## [1,]     1          1          0          0
modelloRWdummy$H
## , , 1
## 
##      [,1]
## [1,]    0
modelloRWdummy$Q
## , , 1
## 
##      [,1] [,2]
## [1,]   10    0
## [2,]    0    1

RW + stagionalità trigonometrica stocastica

Con le dummy vengono delle matrici impegnative perche lunghe quanto il numero dell dummy meno uno piu gli altri elementi del modello, quindi faremo solo l’esempio con 6 dummy stagionali e un random wlak cosi da semplificare la comprensione della componente stagionale nei modelli state space.

\[T=\begin{bmatrix} 1 &0&0 &0&0&0 \\0&\cos\bigg(\frac{2\pi}{periodo} \bigg)& \sin\bigg(\frac{2\pi}{periodo} \bigg) & 0&0&0 \\0&\sin\bigg(\frac{2\pi}{periodo} \bigg) & \cos\bigg(\frac{2\pi}{periodo} \bigg) & 0&0&0\\0&0&0&\cos\bigg(\frac{4\pi}{periodo} \bigg)& \sin\bigg(\frac{4\pi}{periodo} \bigg) & 0 \\0&0&0&\sin\bigg(\frac{4\pi}{periodo} \bigg) & \cos\bigg(\frac{4\pi}{periodo} \bigg) & 0\\ 0&0&0&0&0&\cos\bigg(\frac{6\pi}{periodo} \bigg)\end{bmatrix} ,\\ \space \space Q=\begin{bmatrix} \sigma_\eta^2&0&0&0&0&0\\0&\sigma_\gamma^2&0&0&0&0 \\0&0&\sigma_\gamma^2&0&0&0\\0&0&0&\sigma_\gamma^2&0&0\\0&0&0&0&\sigma_\gamma^2&0\\0&0&0&0&0&\sigma_\gamma^2 \end{bmatrix} , \space \space H=0, \space \space Z=\big[\begin{array}{cc} 1&1&0&1&0 \end{array} \big] \]

N.B. La matrice \(T\) contiene al suo interno la matrice di rotazione:\(R(y)\) Quindi il modello state space diventa:

  • Equazioni di misurazione (o di osservazione):

\[ Y_t = \big[\begin{array}{cc} 1&1&0&1&0\end{array}\big] \begin{bmatrix} \alpha_{t+1} \\ sinous1_{t+1}\\ sinous2_{t+1} \\ sinous3_{t+1}\\ sinous4_{t+1} \\ sinous5_{t+1} \end{bmatrix} \] - Equazioni di stato (o di transizione):

\[ \begin{bmatrix} \alpha_{t+1} \\ sinous1_{t+1}\\ sinous2_{t+1} \\ sinous3_{t+1}\\ sinous4_{t+1} \\ sinous5_{t+1} \end{bmatrix} = \begin{bmatrix} 1 &0&0 &0&0&0 \\0&\cos\bigg(\frac{2\pi}{periodo} \bigg)& \sin\bigg(\frac{2\pi}{periodo} \bigg) & 0&0&0 \\0&\sin\bigg(\frac{2\pi}{periodo} \bigg) & \cos\bigg(\frac{2\pi}{periodo} \bigg) & 0&0&0\\0&0&0&\cos\bigg(\frac{4\pi}{periodo} \bigg)& \sin\bigg(\frac{4\pi}{periodo} \bigg) & 0 \\0&0&0&\sin\bigg(\frac{4\pi}{periodo} \bigg) & \cos\bigg(\frac{4\pi}{periodo} \bigg) & 0\\ 0&0&0&0&0&\cos\bigg(\frac{6\pi}{periodo} \bigg)\end{bmatrix} \begin{bmatrix} \alpha_{t} \\ sinous1_{t}\\ sinous2_{t} \\ sinous3_{t} \\ sinous4_{t} \\ sinous5_{t} \end{bmatrix} +\begin{bmatrix} \nu_t\\\gamma_t \\\gamma_t\\\gamma_t\\\gamma_t\\\gamma_t \end{bmatrix} \space, \\ \space \space con \space \space \nu_t \sim WN(0, \sigma_\eta^2) \\ \space \space \space con \space \space \gamma_t \sim WN(0, \sigma_\gamma^2) \]

library(KFAS)
y<-AirPassengers
sigma2_eta=10
sigma2_epsilon=2
sigma2_k=0.1
sigma2_gamma=1

serie_di_Furier<-2*pi*(1:3)/6
f1<-cos(serie_di_Furier)
f2<-sin(serie_di_Furier)

periodo=6
modelloRWsinus<-SSModel(y~0+SSMtrend(degree = 1,Q=sigma2_eta)+SSMseasonal(period =periodo,Q= sigma2_gamma, sea.type = "trigonometric" ), H= 0)

print("coseno");f1;print("seno");f2
## [1] "coseno"
## [1]  0.5 -0.5 -1.0
## [1] "seno"
## [1] 8.660254e-01 8.660254e-01 1.224606e-16
modelloRWsinus$T
## , , 1
## 
##            level  sea_trig1 sea_trig*1  sea_trig2 sea_trig*2 sea_trig3
## level          1  0.0000000  0.0000000  0.0000000  0.0000000         0
## sea_trig1      0  0.5000000  0.8660254  0.0000000  0.0000000         0
## sea_trig*1     0 -0.8660254  0.5000000  0.0000000  0.0000000         0
## sea_trig2      0  0.0000000  0.0000000 -0.5000000  0.8660254         0
## sea_trig*2     0  0.0000000  0.0000000 -0.8660254 -0.5000000         0
## sea_trig3      0  0.0000000  0.0000000  0.0000000  0.0000000        -1
modelloRWsinus$Z
## , , 1
## 
##      level sea_trig1 sea_trig*1 sea_trig2 sea_trig*2 sea_trig3
## [1,]     1         1          0         1          0         1
modelloRWsinus$H
## , , 1
## 
##      [,1]
## [1,]    0
modelloRWsinus$Q
## , , 1
## 
##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]   10    0    0    0    0    0
## [2,]    0    1    0    0    0    0
## [3,]    0    0    1    0    0    0
## [4,]    0    0    0    1    0    0
## [5,]    0    0    0    0    1    0
## [6,]    0    0    0    0    0    1