UCM Models in State Space Form in KFAS

In-depth Articles

State Space Models in KFAS

General Definition

The state space form used by KFAS is slightly different from what we have used so far. In particular, it is not possible to use the vectors ct and dt, and the covariance matrix of the state vector disturbances is specified using the product of two different matrices. KFAS uses the following state space form:

  • Measurement (or Observation) Equations: establish a linear relationship between a time series (or a vector of time series) and \(m\) state variables that contain the unobservable components in the case of UCMs and other components that help construct them.\[ \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) \] where \(\alpha_t\) is a state vector that does not necessarily have to have the same dimension as the time series.
  • State (or Transition) Equations: inform about the law of motion, that is, how \(\alpha_t\) evolves over time (it is a VAR(1) with coefficients that can evolve over time).\[ \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) \]This form is called the future form because on the left there is only time \(t+1\) and on the right time \(t\). Usually one writes \(\ni_{t+1}\) but since it is a WN it doesn't change anything. The task of matrix \(R\) is to assign the two disturbances to the correct transition equations: for example, associating the first disturbance with the second equation and the second disturbance with the third equation.

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.

Technical conditions:

  • First condition: \(E(\underline\epsilon_t \space \underline\nu_t^T)=G_t\) and it is assumed that \(G_t=0\) for simplicity;
  • Second condition: \(E(\underline\epsilon_t \space \underline\nu_s^T)=0\) for every \(t \neq s\);
  • Third condition: \(E[(\underline\alpha_1-\underline\alpha_{1|0})\space \underline\epsilon_t^T]=0\);
  • Fourth condition: \(E[(\underline\alpha_1-\underline\alpha_{1|0})\space \underline\nu_t^T]=0\).

UCM Definition in KFAS

The R package that we will use for estimating models in state space form is KFAS, which you can download from CRAN like all major R packages. The package contains excellent algorithms for performing inference on state variables (Kalman filter and smoothers). The KFAS package defines a class of objects called SSModel, which is in fact a list containing all the information of the state space form. The constructor, that is the function that generates objects of type SSModel, is itself called SSModel() and, for our purposes, requires only three formal parameters: formula, data, and H. Therefore, to generate a state space form and assign it to the model variable we will write:

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 \] Therefore the state space model becomes:

  • Measurement (or Observation) Equations:

\[ Y_t = \alpha_t + \epsilon_t \space, \space \space con \space \space \epsilon_t \sim WN(0, \sigma_\epsilon^2) \] - State (or Transition) Equations:

\[ \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] \] Therefore the state space model becomes:

  • Measurement (or Observation) Equations:

\[ 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) \] - State (or Transition) Equations:

\[ \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] \] Therefore the state space model becomes:

  • Measurement (or Observation) Equations:

\[ 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} \] - State (or Transition) Equations:

\[ \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)\) Therefore the state space model becomes:

  • Measurement (or Observation) Equations:

\[ 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} \] - State (or Transition) Equations:

\[ \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)\) Therefore the state space model becomes:

  • Measurement (or Observation) Equations:

\[ 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} \] - State (or Transition) Equations:

\[ \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)\) Therefore the state space model becomes:

  • Measurement (or Observation) Equations:

\[ 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} \] - State (or Transition) Equations:

\[ \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