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 matrixs. 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 is completata from the specificazione of the 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 components stazionarie it inseriscono the mean e the variance di tali components, nel caso invece di components not stazionarie (come trend o the seasonaltà) it possono inserire condizioni diffuse, i.e. a mean arbitraria e of the varianze infinite (c’è infinita incertezza with respect al value che the v.c. può assumere). La possibilità di porre condizioni diffuse garantisce che definire \(\alpha_{1|0}\) e \(P_{1|0}\) not sia a limite ma also l’algoritmo usato not the consente it può fissare a variance 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 eestimateting 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:

model <- SSModel(formula, data, H)

Il according to e the terzo parameter formale are very banali: data deve contenere a data.frame con the dati necessarialmodel,cioè the o the variables dependent e the eventuali regressori(se the variables are disponibili nel global environment the parameter formale data può essere omesso); H deve contenere the matrix di covariance (o l’array di matrici di covariance if these evolvono nel tempo) dell’error di observation. Al parameter formale formula deve essere passato a oggetto di tipo formula comune in R per gestire variables esogene ed endogene. Per generare a forma state space partendo direttamente from the matrici di sistema it deve utilizzare all’interno del parameter formula the function:

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

In KFAS esistono diverse functions utilizzabili additivamente per creare the components more comuni of the models a components not osservabili. Per the dettagli of the functions it consiglia di riferirsi to the documentazione del pacchetto KFAS, qui ci limitiamo a mostrare of the esempi per l’utilizzo in model univariati. SSMarima(ar, ma, d = 0, Q) produce the forma state space di a model ARIMA(p,d,q) con the coefficients of the parte AR e MA specificati in the vettori ar e ma, rispettivamente.

Nella prossima parte vedremo nel dettaglio the specificazione of the varie components of the models 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

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


modelRW$T
## , , 1
## 
##       level
## level     1
modelRW$Z
## , , 1
## 
##      level
## [1,]     1
modelRW$H
## , , 1
## 
##      [,1]
## [1,]    2
modelRW$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

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


modelRWD$T
## , , 1
## 
##       level slope
## level     1     1
## slope     0     1
modelRWD$Z
## , , 1
## 
##      level slope
## [1,]     1     0
modelRWD$H
## , , 1
## 
##      [,1]
## [1,]    1
modelRWD$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
modelRWcycle<-SSModel(y~0+SSMtrend(degree = 1,Q=sigma2_eta)+SSMcycle(period =periodo,Q= sigma2_k ), H= 0)

modelRWcycle$T
## , , 1
## 
##        level      cycle    cycle*
## level      1  0.0000000 0.0000000
## cycle      0  0.8412535 0.5406408
## cycle*     0 -0.5406408 0.8412535
modelRWcycle$Z
## , , 1
## 
##      level cycle cycle*
## [1,]     1     1      0
modelRWcycle$H
## , , 1
## 
##      [,1]
## [1,]    0
modelRWcycle$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 matrix \(T\) contiene al suo interno the matrix 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
modelRWDcycle<-SSModel(y~0+SSMtrend(degree = 2,Q=list(sigma2_eta,0))+SSMcycle(period =periodo,Q= sigma2_k ), H= 0)

modelRWDcycle$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
modelRWDcycle$Z
## , , 1
## 
##      level slope cycle cycle*
## [1,]     1     0     1      0
modelRWDcycle$H
## , , 1
## 
##      [,1]
## [1,]    0
modelRWDcycle$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 + seasonaltà a dummy

Con the dummy are of the matrici impegnative perche lunghe quanto the numero dell dummy less uno piu the altri elementi del model, therefore faremo only l’esempio con 4 dummy seasonal e a random wlak cosi da semplificare the comprensione of the component seasonal in the models 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 matrix \(T\) contiene al suo interno the matrix 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
modelRWdummy<-SSModel(y~0+SSMtrend(degree = 1,Q=sigma2_eta)+SSMseasonal(period =periodo,Q= sigma2_gamma, sea.type = "dummy" ), H= 0)

modelRWdummy$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
modelRWdummy$Z
## , , 1
## 
##      level sea_dummy1 sea_dummy2 sea_dummy3
## [1,]     1          1          0          0
modelRWdummy$H
## , , 1
## 
##      [,1]
## [1,]    0
modelRWdummy$Q
## , , 1
## 
##      [,1] [,2]
## [1,]   10    0
## [2,]    0    1

RW + seasonaltà trigonometrica stocastica

Con the dummy are of the matrici impegnative perche lunghe quanto the numero dell dummy less uno piu the altri elementi del model, therefore faremo only l’esempio con 6 dummy seasonal e a random wlak cosi da semplificare the comprensione of the component seasonal in the models 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 matrix \(T\) contiene al suo interno the matrix 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
modelRWsinus<-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
modelRWsinus$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
modelRWsinus$Z
## , , 1
## 
##      level sea_trig1 sea_trig*1 sea_trig2 sea_trig*2 sea_trig3
## [1,]     1         1          0         1          0         1
modelRWsinus$H
## , , 1
## 
##      [,1]
## [1,]    0
modelRWsinus$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