# TSSL Lab 2 - Structural model, Kalman filtering and EM
We will continue to work with the Global Mean Sea Level (GMSL) data that we got acquainted with in lab 1. The data is taken from https://climate.nasa.gov/vital-signs/sea-level/ and is available on LISAM in the file `sealevel.csv`.

In this lab we will analyse this data using a structural time series model. We will first set up a model and implement a Kalman filter to infer the latet states of the model, as well doing long-term prediction. We will then implement a disturbance smoother and an expectation maximization algorithm to tune the parameters of the model. 

We load a few packages that are useful for solving this lab assignment.


```python
import pandas  # Loading data / handling data frames
import numpy as np
import scipy
import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (12,8)  # Increase default size of plots
```

## 2.1 Setting up a structural state space model

We start by loading and plotting data to remind ourselves what it looks like.


```python
data=pandas.read_csv('sealevel.csv',header=0)
```


```python
y = data['GMSL'].values
u = data['Year'].values
ndata = len(y)

plt.plot(figsize=(12,8))
plt.plot(u,y)
plt.xlabel('Year')
plt.ylabel('GMSL deviation (mm)')
plt.show()
```

In this lab we will use a structural time series model to analys this data set. Specifically, we assume that the data $\{y_t\}_{t\geq 1}$ is generated by

\begin{align}
    y_t = \mu_t + \gamma_t + \varepsilon_t
\end{align}

where $\mu_t$ is a trend component, $\gamma_t$ is a seasonal component, and $\varepsilon_t$ is an observation noise. The model is expressed using a state space representation,

\begin{align}
    \alpha_{t+1} &= T \alpha_t + R\eta_t, & \eta_t&\sim N(0,Q), \\
    y_t &= Z \alpha_t + \varepsilon_t, & \varepsilon_t&\sim N(0,\sigma_\varepsilon^2).
\end{align}

**Q0:** Let $d = \dim(\alpha_t)$ denote the _state dimension_ and $d_\eta = \dim(\eta_t)$ denote the dimension of the state noise. Then, what are the dimenisons of the matrices $T$, $R$, and $Z$ of the state space model?

**A:** dim(T) = (k+s-1)X(k+s-1) , dim(R) = (k+s-1)x2 , dim(Z) = 1x(k+s-1)

**Q1:** Create the state space matrices $T_{[\mu]}$, $R_{[\mu]}$, and $Z_{[\mu]}$ corresponding to the trend component $\mu_t$. We assume a local linear trend (that is, of order $k=2$). 

_Hint:_ Use **2-dimensional** `numpy.ndarray`s of the correct sizes to represent all the matrices.


```python
k=2
import math as m
T_mu = np.zeros(shape=(k,k))
R_mu = np.zeros(shape=(k,1))
Z_mu = np.zeros(shape=(1,k))
for i in range(k):
    if i == 0:
        for j in range(k):
            t = m.factorial(k)/(m.factorial(k-j-1)*m.factorial(j+1))
            #print("\n",j,m.factorial(k),m.factorial(k-j-1),m.factorial(j+1),t)
            T_mu[i,j] = (-1)**(j+2)*t
    else:
        T_mu[i,i-1] = 1
R_mu[0] = 1
Z_mu[0] = 1
```

**Q2:** There is a yearly seasonal pattern present in the data. What should we set the periodicity $s$ of the seasonal component to, to capture this pattern?

_Hint:_ Count the average number of observations per (whole) year and round to the closest integer.


```python
years = u[:-2].astype(int)
years = np.unique(years)
s = round(len(u)/len(years))
print("The period of the seasonality component is",s)
```

    The period of the seasonality component is 37


**Q3:** What is the _state dimension_ of a seasonal component with periodicity $s$? That is, how many states are needed in the corresponding state space representation?

**A:** dim(T_gamma)=(s-1)x(s-1) dim(R_gamma) = (s-1)x1 dim(Z_gamma) = 1x(s-1)


**Q4:** Create the state space matrices $T_{[\gamma]}$, $R_{[\gamma]}$, and $Z_{[\gamma]}$ corresponding to the seasonal component $\gamma_t$. 

_Hint:_ Use **2-dimensional** `numpy.ndarray`s of the correct sizes to represent all the matrices.


```python
T_gamma = np.zeros(shape=(s-1,s-1))
R_gamma = np.zeros(shape=(s-1,1))
Z_gamma = np.zeros(shape=(1,s-1))
c = T_gamma.shape[1]
for i in range(c):
    if(i==0):
        T_gamma[i,]=-1
    else:
        T_gamma[i,i-1] = 1
R_gamma[0,] = 1
Z_gamma[0,0] = 1
#print("Tgamma",T_gamma)
#print("Rgamma",R_gamma)
#print("Zgamma",Z_gamma)
```

**Q5:** Using the matrices that you have constructed above, create the state space matrices for the complete structural time series model. Print out the shapes of the resulting system matrices and check that they correspond to what you expect (cf **Q0**).

_Hint:_ Use `scipy.linalg.block_diag` and `numpy.concatenate`.


```python
T = np.concatenate((np.concatenate((T_mu,np.zeros(shape=(k,s-1))),axis=1),
                    np.concatenate((np.zeros(shape=(s-1,k)),T_gamma),axis=1)),axis=0)
R = np.concatenate((np.column_stack((R_mu,np.zeros(shape=(k,1)))),
                    np.column_stack((np.zeros(shape=(s-1,1)),R_gamma))),axis=0)
Z = np.append(Z_mu,Z_gamma,axis=1)
#print("T",T)
#print("R",R)
#print("Z",Z)
```

We also need to specify the variances of the process noise $\eta_t$ and measurement noise $\varepsilon_t$. Below, we will estimate (two of) these variances from data, but for now we set them arbitrarily to get an initial model to work with.


```python
# Some arbitrary noise values for now
sigma_trend = 0.01
sigma_seas = 1
sigma_eps = 1
Q = np.array([[sigma_trend**2, 0.], [0., sigma_seas**2]])  # Process noise covariance matrix
#print(Q)
```

Finally, to complete the model we need to specify the distribution of the initial state. This encodes our _a priori_ belief about the actual values of the trend and seasonality, i.e., before observing any data.

**Q6:** Set up the mean vector of the initial state $a_1 = \mathbb{E}[\alpha_1]$ such that:
* The trend component starts at the first observation, $\mathbb{E}[\mu_1] = y_1$,
* The slope of the trend is _a priori_ zero in expectation, $\mathbb{E}[\mu_1 - \mu_0] = 0$,
* The initial mean of all states related to the seasonal component are zero.

Also, create an initial state covariance matrix $P_1 = \text{Cov}(\alpha_1)$ as an identity matrix of the correct dimension, multiplied with a large value (say, 100) to represent our uncertainty about the initial state.


```python
a1 = np.zeros(shape=(s+k-1,1))
a1[0] = y[0]
a1[1] = y[0]
P1 = 100*np.identity(s+k-1)
#print("alpha1",alpha1)
#print("P1",P1)
```

We have now defined all the matrices etc. that make up the structural state space model. For convenience, we can create an object of the class `LGSS` available in the module `tssltools_lab2` as a container for these quantities.


```python
from tssltools_lab2 import LGSS  # Module available in LISAM
model = LGSS(T, R, Q, Z, sigma_eps**2, a1, P1)

help(model.get_params)
```

    Help on method get_params in module tssltools_lab2:
    
    get_params() method of tssltools_lab2.LGSS instance
        Return all model parameters.
        
        T, R, Q, Z, H, a1, P1 = model.get_params()
    


## 2.2 Kalman filtering for the structural model
Now we have the data and a model available. Next, we will turn our attention to the inference problem, which is a central task when analysing time series data using the state space framework.

State inference is the problem of estimating the unknown (latent) state variables given the data. For the time being we assume that the _model parameters_ are completely specified, according to above, and only consider how to estimate the states using the Kalman filter.

In the questions below we will treat the first $n=800$ time steps as training data and the remaining $m=197$ observations as validation data. 


```python
n = 800
m = ndata-n
```


**Q7:** Complete the Kalman filter implementation below. The function should be able to handle missing observations, which are encoded as "not a number", i.e. `y[t] = np.nan` for certain time steps `t`. 

_Hint:_ The Kalman filter involves a lot of matrix-matrix and matrix-vector multiplications. It turns out to be convient to store sequences of vectors (such as the predicted and filtered state estimates) as `(d,1,n)` arrays, instead of `(d,n)` or `(n,d)` arrays. In this way the matrix multiplications will result in 2d-arrays of the correct shapes without having to use a lot of explicit `reshape`. However, clearly, this is just a matter of coding style preferences!


```python
from tssltools_lab2 import kfs_res

def kalman_filter(y, model: LGSS):
    """Kalman filter for LGSS model with one-dimensional observation.

    :param y: (n,) array of observations. May contain nan, which encodes missing observations.
    :param model: LGSS object with the model specification.
    
    :return kfs_res: Container class with member variables,
        alpha_pred: (d,1,n) array of predicted state means.
        P_pred: (d,d,n) array of predicted state covariances.
        alpha_filt: (d,1,n) array of filtered state means.
        P_filt: (d,d,n) array of filtered state covariances.
        y_pred: (n,) array of means of p(y_t | y_{1:t-1})
        F_pred: (n,) array of variances of p(y_t | y_{1:t-1})
    """

    n = len(y)
    d = model.d  # State dimension
    alpha_pred = np.zeros((d, 1, n))
    P_pred = np.zeros((d, d, n))
    alpha_filt = np.zeros((d, 1, n))
    P_filt = np.zeros((d, d, n))
    y_pred = np.zeros(n)
    F_pred = np.zeros(n)

    T, R, Q, Z, H, a1, P1 = model.get_params()  # Get all model parameters (for brevity)

    for t in range(n):
        # Time update (predict)
        if t == 0:  # Initialize predictions at first time step
            alpha_pred[:,:,t] = a1
            P_pred[:,:,t] = P1
        else:  # All consecutive time steps
            alpha_pred[:,:,t] = T @ alpha_filt[:,:,t-1]
            P_pred[:,:,t] = T @ P_filt[:,:,t-1] @ T.T + R @ Q @ R.T

        # Compute prediction of current output
        y_pred[t] = Z @ alpha_pred[:,:,t] # complete this line
        F_pred[t] = Z @ P_pred[:,:,t] @ Z.T + H  # complete this line

        # Measurement update (correct)
        if np.isnan(y[t]): # Handle missing data
            alpha_filt[:,:,t] = alpha_pred[:,:,t]
            P_filt[:,:,t] = P_pred[:,:,t]
        else:
            k_t = P_pred[:,:,t] @ Z.T / F_pred[t]
            alpha_filt[:,:,t] = alpha_pred[:,:,t] + (k_t * (y[t] - y_pred[t]))
            P_filt[:,:,t] = (np.identity(d) - k_t @ Z) @ P_pred[:,:,t]

    kf = kfs_res(alpha_pred, P_pred, alpha_filt, P_filt, y_pred, F_pred)
    return kf
```

**Q8:** Use the Kalman filter to infer the states of the structural time series applied to the sealevel data. Run the filter on the training data (i.e., first $n=700$ time steps), followed by a long-range prediction of $y_t$ for the remaining time points. 

Generate a plot which shows:
1. The data $y_{1:n+m}$,
2. The one-step predictions $\hat y_{t|t-1} \pm 1 $ standard deviation for the training data, i.e., $t = 1,...,n$,
3. The long-range predictions $\hat y_{t|n} \pm 1 $ standard deviation for the validation data, i.e., $t= n+1,...,n+m$,
4. A vertical line indicating the switch between training and validation data, using, e.g., `plt.axvline(x=u[n])`.

_Hint:_ It is enough to call the `kalman_filter` function once. Make use of the missing data functionality!


```python
a = np.empty(m)
a[:] = np.NaN
ynew = np.append(y[0:800],a)
kf_res = kalman_filter(ynew,model)
y_pred = kf_res.y_pred
errors = kf_res.F_pred
upper = y_pred + np.sqrt(errors)
lower = y_pred - np.sqrt(errors)
```


```python
plt.plot(u,y,label="Actual Data")
plt.plot(u,y_pred,c='yellow',label="Predictions")
plt.plot(u,upper,c='pink',label="Upper CI")
plt.plot(u,lower,c='c',label="Lower CI")
plt.axvline(x=u[n],c='r',label="Switch")
plt.legend()
plt.xlabel('Year')
plt.ylabel('GMSL deviation (mm)')
plt.show()
```

**Q9:** Based on the output of the Kalman filter, compute the training data log-likelihood $\log p(y_{1:n})$.


```python
def llik(y,y_pred,errors) :
    n = len(y)
    return(-(n/2)*np.log(2*np.pi)-(1/2)*sum(np.log(errors)+((y-y_pred)**2)/errors))
print("The log-likelyhood is ",llik(y[0:(n-1)],y_pred[0:(n-1)],errors[0:(n-1)]))
```

    The log-likelyhood is  -2626.3783139837383


## 2.3 Identifying the noise variances using the EM algorithm
So far we have used fixed model parameters when running the filter. In this section we will see how the model parameters can be learnt from data using the EM algorithm. Specifically, we will try to learn the variance of the state noise affecting the seasonal component as well as the variance of the observation noise,

\begin{align}
    \theta = (\sigma_\gamma^2, \sigma_\varepsilon^2).
\end{align}

For brevity, the variance of the trend component $\sigma_\mu^2$ is fixed to the value $\sigma_\mu^2 = 0.01^2$ as above. (See Appendix A below for an explanation.)

Recall that we consider $y_{1:n}$ as the training data, i.e., we will estimate $\theta$ using only the first $n=800$ observations.

**Q10:** Which optimization problem is it that the EM algorithm is designed to solve? Complete the line below!

**A:** $\hat\theta = \arg\max_{\theta} Q(\theta,\tilde \theta) = const-\frac{1}{2}\sum\limits_{t=1}^{n}[log|\sigma_{\epsilon}^2|+log|Q|+(\hat\epsilon_{t|n}^2+var(\epsilon_{t}|y_{1:n}))\sigma_{\epsilon}^{-2}+trace[(\hat\eta_{t|n}{\hat\eta_{t|n}}^T+var(\eta|y_{1:n}))Q^{-1}] $

**Q11:** Write down the updating equations on closed form for the M-step in the EM algorithm.

**A:** ${\hat\sigma_{\epsilon}}^2 = \frac{1}{n} \sum\limits_{t=1}^{n}(\hat\epsilon_{t|n}^2+var(\epsilon_t|y_{1:n}))$
$\hat\sigma_\gamma^2 =  \frac{1}{n} \sum\limits_{t=1}^{n}d_t/n$ 

$d_t=((\hat\eta_{t|n}{\hat\eta_{t|n}}^T+var(\eta|y_{1:n})))[n,n]$

To implement the EM algorithm we need to solve a _smoothing problem_. The Kalman filter that we implemented above is based only on a forward propagation of information. The _smoother_ complements the forward filter with a backward pass to compute refined state estimates. Specifically, the smoothed state estimates comprise the mean and covariances of

\begin{align}
    &p(\alpha_t \mid y_{1:n}), & t=&1,\dots,n
\end{align}

Furthermore, when implementing the EM algorithm it is convenient to work with the (closely related) smoothed estimates of the disturbances, i.e., the state and measurement noise,

\begin{align}
    &p(\eta_t \mid y_{1:n}), & t=&1,\dots,n-1 \\
    &p(\varepsilon_t \mid y_{1:n}), & t=&1,\dots,n
\end{align}

An implementation of a state and disturbance smoother is available in the `tssltools_lab2` module. You may use this when implementing the EM algorithm below.



```python
from tssltools_lab2 import kalman_smoother
help(kalman_smoother)
```

    Help on function kalman_smoother in module tssltools_lab2:
    
    kalman_smoother(y, model: tssltools_lab2.LGSS, kf: tssltools_lab2.kfs_res)
        Kalman (state and disturbance) smoother for LGSS model with one-dimensional observation.
        
        :param y: (n,) array of observations. May contain nan, which encodes missing observations.
        :param model: LGSS object with the model specification.
        :parma kf: kfs_res object with result from a Kalman filter foward pass.
        
        :return kfs_res: Container class. The original Kalman filter result is augmented with the following member variables,
            alpha_sm: (d,1,n) array of smoothed state means.
            V: (d,d,n) array of smoothed state covariances.
            eps_hat: (n,) array of smoothed means of observation disturbances.
            eps_var: (n,) array of smoothed variances of observation disturbances.
            eta_hat: (deta,1,n) array of smoothed means of state disturbances.
            eta_cov: (deta,deta,n) array of smoothed covariances of state disturbances.
    


**Q12:** Implement an EM algorithm by completing the code below. Run the algorithm for 100 iterations and plot the traces of the parameter estimates, i.e., the values $\theta_r$, for $r = 0,\dots,100$.

_Note:_ When running the Kalman filter as part of the EM loop you should only filter the _training data_ (i.e. excluding the prediction for validation data).



```python
num_iter = 100

def compute_theta(ks):
    """Implements the M-step of the EM-algorithm, based on the Kalman smoother results (ks). See Q11."""
    eta_hat = ks.eta_hat
    eta_cov = ks.eta_cov
    eps_hat = ks.eps_hat
    eps_var = ks.eps_var
    
    n = len(eps_hat)
    
    sigma2_eps_hat = (1/n)*sum(eps_hat**2+eps_var)
    
    sigma2_gamma = 0 
    for t in range(n):
        d_t = ((eta_hat[:,:,t] @ eta_hat[:,:,t].T)+eta_cov[:,:,t])
        d_t = d_t[1,1]
        sigma2_gamma = sigma2_gamma + d_t
    sigma2_gamma=sigma2_gamma/n
    
    theta_new = [np.sqrt(sigma2_gamma),np.sqrt(sigma2_eps_hat)]
    return theta_new

def EM(y,num_iter):
    sigma_trend = 0.01
    sigma_seas = 1
    sigma_eps = 1

    like = np.zeros(shape=(num_iter,1))
    
    trace = np.zeros(shape=(num_iter,2))
    
    for r in range(num_iter):
        # E-Step

        # updating the model parameters
        Q = np.array([[sigma_trend**2, 0.], [0., sigma_seas**2]])  # Process noise covariance matrix
        model = LGSS(T, R, Q, Z, sigma_eps**2, a1, P1)

        # Calculating Kalman filter and Kalman Smoother
        kf_res = kalman_filter(y, model)
        ks = kalman_smoother(y,model,kf_res)

        # Calculating the likelihood
        y_pred = kf_res.y_pred
        errors = kf_res.F_pred
        like[r,0] = llik(y,y_pred,errors)
        if r!=0:
            if abs(like[r]-like[r-1])<0.01:
                break
        # M-step
        theta = compute_theta(ks)
        sigma_seas = theta[0]
        sigma_eps = theta[1]
        
        trace[r,0]=theta[0]
        trace[r,1]=theta[1]

        print("The theta in iteration is ",r,theta,"with likelihood",like[r,0])
        
    EM_res = [model,like,kf_res,trace]
        
    return(EM_res)
em_res=EM(y,100)
```

    The theta in iteration is  0 [1.1756464079272044, 1.5587868045043296] with likelihood -3221.954007852065
    The theta in iteration is  1 [1.154401913561508, 1.9376390158770391] with likelihood -2775.4012571485014
    The theta in iteration is  2 [1.0919471486133427, 2.1681479176554164] with likelihood -2688.059990725837
    The theta in iteration is  3 [1.027022512719855, 2.293804013068092] with likelihood -2661.1132721657473
    The theta in iteration is  4 [0.967825435006054, 2.358840340125723] with likelihood -2649.191987032433
    The theta in iteration is  5 [0.915348291082476, 2.3931392029523577] with likelihood -2641.5461853192096
    The theta in iteration is  6 [0.8689750461211504, 2.412908443927547] with likelihood -2635.541353613631
    The theta in iteration is  7 [0.8278225249690394, 2.425955268879042] with likelihood -2630.4745310876965
    The theta in iteration is  8 [0.791069043266986, 2.4357992949840823] with likelihood -2626.088732026287
    The theta in iteration is  9 [0.7580225227018894, 2.4439760206737007] with likelihood -2622.242436266194
    The theta in iteration is  10 [0.7281176152261115, 2.4511549223632847] with likelihood -2618.8368973616984
    The theta in iteration is  11 [0.7008955994974848, 2.4576403967841953] with likelihood -2615.7971153669114
    The theta in iteration is  12 [0.6759831538061377, 2.463584877142113] with likelihood -2613.0643569578915
    The theta in iteration is  13 [0.6530743390567524, 2.469076439323705] with likelihood -2610.5918753113206
    The theta in iteration is  14 [0.6319163500454358, 2.4741741802123323] with likelihood -2608.342005816739
    The theta in iteration is  15 [0.6122985138519278, 2.47892257533545] with likelihood -2606.284077128845
    The theta in iteration is  16 [0.594043834762242, 2.483357543551083] with likelihood -2604.3928713187324
    The theta in iteration is  17 [0.5770024710289011, 2.4875092449817258] with likelihood -2602.6474707803263
    The theta in iteration is  18 [0.5610466643802021, 2.4914035561031285] with likelihood -2601.030383507352
    The theta in iteration is  19 [0.5460667642773015, 2.4950629700931355] with likelihood -2599.526871836943
    The theta in iteration is  20 [0.5319680827593796, 2.498507216250323] with likelihood -2598.1244318740714
    The theta in iteration is  21 [0.5186683848928374, 2.501753720852812] with likelihood -2596.812385832498
    The theta in iteration is  22 [0.506095870022656, 2.5048179660828054] with likelihood -2595.581559914157
    The theta in iteration is  23 [0.49418753538432486, 2.5077137771026963] with likelihood -2594.4240276392356
    The theta in iteration is  24 [0.48288784012363306, 2.5104535555725187] with likelihood -2593.332903724343
    The theta in iteration is  25 [0.4721476072276734, 2.5130484719116457] with likelihood -2592.3021773391993
    The theta in iteration is  26 [0.4619231152914508, 2.5155086251280188] with likelihood -2591.326576290055
    The theta in iteration is  27 [0.45217534283187727, 2.517843176786272] with likelihood -2590.401455677614
    The theta in iteration is  28 [0.442869336002117, 2.5200604641242883] with likelihood -2589.522706062277
    The theta in iteration is  29 [0.43397367675539184, 2.5221680961982322] with likelihood -2588.6866772828894
    The theta in iteration is  30 [0.42546003326099, 2.5241730360981967] with likelihood -2587.890114916715
    The theta in iteration is  31 [0.41730277804974863, 2.526081671644245] with likelihood -2587.1301070098016
    The theta in iteration is  32 [0.40947866222699747, 2.5278998764892657] with likelihood -2586.4040391992226
    The theta in iteration is  33 [0.40196653633316376, 2.529633063181801] with likelihood -2585.70955672972
    The theta in iteration is  34 [0.394747110200926, 2.5312862294501293] with likelihood -2585.0445321634825
    The theta in iteration is  35 [0.3878027455614745, 2.5328639987396824] with likelihood -2584.407037814503
    The theta in iteration is  36 [0.38111727627286823, 2.5343706558541537] with likelihood -2583.795322121906
    The theta in iteration is  37 [0.3746758519427221, 2.535810178404161] with likelihood -2583.207789321886
    The theta in iteration is  38 [0.3684648014430029, 2.537186264650858] with likelihood -2582.642981894049
    The theta in iteration is  39 [0.3624715134030326, 2.5385023582364257] with likelihood -2582.0995653503796
    The theta in iteration is  40 [0.3566843312461034, 2.539761670215544] with likelihood -2581.576315010335
    The theta in iteration is  41 [0.35109246072744543, 2.5409671987385387] with likelihood -2581.0721044657766
    The theta in iteration is  42 [0.3456858882537856, 2.5421217466841823] with likelihood -2580.5858954890755
    The theta in iteration is  43 [0.3404553085309849, 2.543227937496143] with likelihood -2580.1167291776383
    The theta in iteration is  44 [0.33539206030691776, 2.544288229441342] with likelihood -2579.6637181615533
    The theta in iteration is  45 [0.3304880691603748, 2.5453049284771385] with likelihood -2579.2260397279474
    The theta in iteration is  46 [0.3257357964401206, 2.5462801998892974] with likelihood -2578.8029297382845
    The theta in iteration is  47 [0.3211281935867786, 2.547216078840473] with likelihood -2578.3936772337047
    The theta in iteration is  48 [0.3166586611783157, 2.5481144799508835] with likelihood -2577.997619638753
    The theta in iteration is  49 [0.31232101213111974, 2.5489772060171743] with likelihood -2577.6141384872776
    The theta in iteration is  50 [0.30810943856589174, 2.5498059559616086] with likelihood -2577.2426556048295
    The theta in iteration is  51 [0.30401848191314984, 2.550602332092987] with likelihood -2576.8826296914785
    The theta in iteration is  52 [0.3000430058889881, 2.5513678467501033] with likelihood -2576.5335532564227
    The theta in iteration is  53 [0.2961781720194421, 2.5521039283907414] with likelihood -2576.194949862632
    The theta in iteration is  54 [0.29241941743266747, 2.5528119271811254] with likelihood -2575.8663716451597
    The theta in iteration is  55 [0.2887624346732283, 2.5534931201348945] with likelihood -2575.5473970716366
    The theta in iteration is  56 [0.28520315332301055, 2.554148715845039] with likelihood -2575.237628917536
    The theta in iteration is  57 [0.28173772323934587, 2.5547798588469837] with likelihood -2574.936692432301
    The theta in iteration is  58 [0.27836249924351, 2.555387633647393] with likelihood -2574.644233675316
    The theta in iteration is  59 [0.2750740271123259, 2.555973068449096] with likelihood -2574.359918003523
    The theta in iteration is  60 [0.271869030742621, 2.5565371385995723] with likelihood -2574.0834286944532
    The theta in iteration is  61 [0.26874440037311054, 2.5570807697872127] with likelihood -2573.8144656905342
    The theta in iteration is  62 [0.26569718176122836, 2.5576048410076178] with likelihood -2573.5527444522277
    The theta in iteration is  63 [0.2627245662237413, 2.5581101873193584] with likelihood -2573.297994908858
    The theta in iteration is  64 [0.25982388145992236, 2.5585976024071106] with likelihood -2573.0499604974316
    The theta in iteration is  65 [0.25699258308477413, 2.5590678409680727] with likelihood -2572.808397280743
    The theta in iteration is  66 [0.2542282468074687, 2.5595216209361524] with likelihood -2572.573073137062
    The theta in iteration is  67 [0.25152856119694206, 2.559959625556944] with likelihood -2572.343767014577
    The theta in iteration is  68 [0.24889132098255523, 2.560382505325434] with likelihood -2572.1202682444527
    The theta in iteration is  69 [0.24631442084302846, 2.5607908797969023] with likelihood -2571.902375907047
    The theta in iteration is  70 [0.24379584964153667, 2.5611853392812343] with likelihood -2571.6898982464745
    The theta in iteration is  71 [0.24133368506902647, 2.56156644642902] with likelihood -2571.482652129033
    The theta in iteration is  72 [0.2389260886615107, 2.561934737717911] with likelihood -2571.280462541689
    The theta in iteration is  73 [0.23657130116039427, 2.562290724846283] with likelihood -2571.083162127021
    The theta in iteration is  74 [0.23426763818782984, 2.562634896041369] with likelihood -2570.890590751538
    The theta in iteration is  75 [0.23201348621172563, 2.56296771728755] with likelihood -2570.702595104372
    The theta in iteration is  76 [0.22980729877738335, 2.563289633480876] with likelihood -2570.519028324004
    The theta in iteration is  77 [0.22764759298484533, 2.563601069514668] with likelihood -2570.3397496504535
    The theta in iteration is  78 [0.22553294619292552, 2.5639024313012064] with likelihood -2570.1646241009976
    The theta in iteration is  79 [0.22346199293259336, 2.564194106733792] with likelihood -2569.993522167461
    The theta in iteration is  80 [0.2214334220139087, 2.5644764665930286] with likelihood -2569.826319533329
    The theta in iteration is  81 [0.21944597381208147, 2.564749865401445] with likelihood -2569.662896809177
    The theta in iteration is  82 [0.2174984377194733, 2.565014642229438] with likelihood -2569.503139284907
    The theta in iteration is  83 [0.21558964975147527, 2.5652711214560493] with likelihood -2569.346936697599
    The theta in iteration is  84 [0.21371849029521373, 2.5655196134872247] with likelihood -2569.1941830137075
    The theta in iteration is  85 [0.21188388199095, 2.565760415434515] with likelihood -2569.044776224594
    The theta in iteration is  86 [0.21008478773687506, 2.5659938117563823] with likelihood -2568.8986181543128
    The theta in iteration is  87 [0.20832020880874993, 2.5662200748649537] with likelihood -2568.7556142789176
    The theta in iteration is  88 [0.20658918308653623, 2.5664394656997436] with likelihood -2568.6156735562818
    The theta in iteration is  89 [0.20489078338077196, 2.566652234270934] with likelihood -2568.4787082658313
    The theta in iteration is  90 [0.2032241158520289, 2.56685862017361] with likelihood -2568.3446338573626
    The theta in iteration is  91 [0.20158831851729495, 2.5670588530750607] with likelihood -2568.2133688084627
    The theta in iteration is  92 [0.19998255983760083, 2.567253153176449] with likelihood -2568.0848344898
    The theta in iteration is  93 [0.19840603738164042, 2.567441731650519] with likelihood -2567.958955037853
    The theta in iteration is  94 [0.1968579765605285, 2.567624791056801] with likelihood -2567.8356572345388
    The theta in iteration is  95 [0.19533762942919877, 2.567802525735409] with likelihood -2567.7148703933053
    The theta in iteration is  96 [0.19384427355027564, 2.5679751221808527] with likelihood -2567.596526251227
    The theta in iteration is  97 [0.19237721091655507, 2.5681427593969186] with likelihood -2567.480558866781
    The theta in iteration is  98 [0.19093576692850858, 2.5683056092336742] with likelihood -2567.366904522896
    The theta in iteration is  99 [0.1895192894234796, 2.5684638367077244] with likelihood -2567.2555016349693



```python
trace = em_res[3]
like = em_res[1]

plt.plot(list(range(100)),trace[:,0],label="sigma_gamma")
plt.plot(list(range(100)),trace[:,1],c='red',label="sigma_eps")
#plt.plot(list(range(100)),like[:,0],c='pink',label="Likelihood")
plt.legend()
plt.xlabel('Year')
plt.ylabel('GMSL deviation (mm)')
plt.show()

plt.plot(list(range(100)),like[:,0],c='m',label="Likelihood")
plt.legend()
plt.xlabel('Year')
plt.ylabel('GMSL deviation (mm)')
plt.show()
```

## 2.4 Further analysing the data
We will now fix the model according to the final output from the EM algorithm and further analyse the data using this model.

**Q13:** Rerun the Kalman filter to compute a _long range prediction for the validation data points,_ analogously to **Q8** (you can copy-paste code from that question). That is, generate a plot which shows:
1. The data $y_{1:n+m}$,
2. The one-step predictions $\hat y_{t|t-1} \pm 1 $ standard deviation for the training data, i.e., $t = 1,...,n$,
3. The long-range predictions $\hat y_{t|n} \pm 1 $ standard deviation for the validation data, i.e., $t= n+1,...,n+m$,
4. A vertical line indicating the switch between training and validation data, using, e.g., `plt.axvline(x=u[n])`.

Furthermore, compute the training data log-likelihood $\log p(y_{1:n})$ using the estimated model (cf. **Q9**).


```python
# updating the model with EM tranined parameters
model = em_res[0]

# Calculating Kalman filter 
kf_res_EM = kalman_filter(ynew,model)
y_pred_EM = kf_res_EM.y_pred
errors_EM = kf_res_EM.F_pred
upper_EM = y_pred + np.sqrt(errors_EM)
lower_EM = y_pred - np.sqrt(errors_EM)

like_filter = llik(y,y_pred_EM,errors_EM)
print("The likelyhood of the filter",like_filter)
```

    The likelyhood of the filter -2803.5722763637373



```python
plt.plot(u,y,label="Actual Data")
plt.plot(u,y_pred_EM,c='yellow',label="Predictions")
plt.plot(u,upper_EM,c='pink',label="Upper CI")
plt.plot(u,lower_EM,c='c',label="Lower CI")
plt.axvline(x=u[n],c='r',label="Switch")
plt.legend()
plt.xlabel('Year')
plt.ylabel('GMSL deviation (mm)')
plt.show()
```

Note that we can view the model for the data $y_t$ as being comprised of an underlying "signal", $s_t = \mu_t + \gamma_t$ plus observation noise $\varepsilon_t$

\begin{align}
    y_t = s_t + \varepsilon_t
\end{align}

We can obtain refined, _smoothed,_ estimates of this signal by conditioning on all the training data $y_{1:n}$. 

**Q14:** Run a Kalman smoother to compute smoothed estimates of the signal, $\mathbb{E}[s_t | y_{1:n}]$, conditionally on all the _training data_. Then, similarly to above, plot the following:
1. The data $y_{1:n+m}$,
2. The smoothed estimates $\mathbb{E}[s_t | y_{1:n}] \pm 1 $ standard deviation for the training data, i.e., $t = 1,...,n$,
3. The predictions $\mathbb{E}[s_t | y_{1:n}] \pm 1 $ standard deviation for the validation data, i.e., $t = n+1,...,n+m$,
4. A vertical line indicating the switch between training and validation data, using, e.g., `plt.axvline(x=u[n])`.

_Hint:_ Express $s_t$ in terms of $\alpha_t$. Based on this expression, compute the smoothed mean and variance of $s_t$ based on the smoothed mean and covariance of $\alpha_t$.


```python
ks_res = kalman_smoother(ynew,model,kf_res_EM)

y_est_pred_sm = np.zeros(len(ynew))
upper = np.zeros(len(ynew))
lower = np.zeros(len(ynew))

for t in range(len(ynew)): 
    y_est_pred_sm[t] = Z @ ks_res.alpha_sm[:,:,t]
    upper[t] = y_est_pred_sm[t] + ks_res.eps_var[t,]
    lower[t] = y_est_pred_sm[t] - ks_res.eps_var[t,]

like_filter = llik(y,y_est_pred_sm,ks_res.eps_var[t,])
print("The likelyhood of the filter",like_filter)
```

    The likelyhood of the filter -4619.272528491351



```python
plt.plot(u,y,label="Actual Data")
plt.plot(u,y_est_pred_sm,c='yellow',label="Predictions")
plt.plot(u,upper,c='pink',label="Upper CI")
plt.plot(u,lower,c='c',label="Lower CI")
plt.axvline(x=u[n],c='r',label="Switch")
plt.legend()
plt.xlabel('Year')
plt.ylabel('GMSL deviation (mm)')
plt.show()
```

**Q15:** Explain, using a few sentences, the qualitative differences (or similarities) between the Kalman filter predictions plotted in **Q13** and the smoothed signal estimates plotted in **Q14** for,
1. Training data points, $t \leq n$
2. Validation data points, $t > n$

**A:** For the traning data points, the smoother knows the behaviour of the data previous and also the future data hence exactly followes even when there is increase/ flip in the peak in each seasonality in the data while the filter does not capture this as it knows only the previous behaviour of the data. Hence the smoother a better choice here.

In the validation data both the smoother and filter do not follow the data in the first seasonality but catch up from the next seasonality. The increase/ flip in data is not observed in validation data, if this was present, the behaviour of smoother and filter can be analyzed in this situaltion which gives a better understanding of which one is better when the data is "unseen" by the algorithm.

We can shed additional light on the properties of the process under study by further decomposing the signal into its trend and seasonal components.

**Q16:** Using the results of the state smoother, compute and plot the _smoothed estimates_ of the two signal components, i.e.:

1. Trend: $\hat \mu_{t|n} = \mathbb{E}[\mu_t | y_{1:n}]$ for $t = 1,\dots,n$
2. Seasonal: $\hat \gamma_{t|n} = \mathbb{E}[\gamma_t | y_{1:n}]$ for $t = 1,\dots,n$

_(You don't have to include confidence intervals here if don't want to, for brevity.)_


```python
mu = np.zeros(shape=(len(ynew),))
gamma = np.zeros(shape=(len(ynew),))
for t in range(1,len(ynew)):
    mu[t] = ks_res.alpha_sm[0,:,t-1]
    gamma[t] = ks_res.alpha_sm[2,:,t-1]
```


```python
plt.plot(u,mu,c='m',label="Mean of Trend Component")
plt.plot(u,gamma,c='yellow',label="Mean of Seasonality Component")
plt.legend()
plt.xlabel('Year')
plt.ylabel('GMSL deviation (mm)')
plt.show()
```

## 2.5 Missing data
We conclude this section by illustrating one of the key merits of the state space approach to time series analysis, namely the simplicity of handling missing data. To this end we will assume that a chunk of observations in the middle of the training data is missing.

**Q17:** Let the values $y_{t}$ for $ 300 < t \leq 400$ be missing (set to `np.nan`). Modify the data and rerun the Kalman filter and smoother. Plot,
1. The Kalman filter predictions, analogously to **Q8**
2. The Kalman smoother predictions, analogously to **Q13**

Comment on the qualitative differences between the filter and smoother estimates and explain what you see (in a couple of sentences).

**A:** The smoother performance has improved and it follows the data better after nans were added to the data, while the filter performance remains the same.


```python
ynan = y.copy()
ynan[300:400] = np.NaN

# Kalman Filter
kf_res_nan = kalman_filter(ynan,model)

y_pred_nan = kf_res_nan.y_pred
errors_nan = kf_res_nan.F_pred
upper_nan = y_pred_nan + np.sqrt(errors)
lower_nan = y_pred_nan - np.sqrt(errors)


#Kalman Smoother
ks_res_nan = kalman_smoother(ynan,model,kf_res_EM)

y_est_pred_sm_nan = np.zeros(len(ynan))
upper = np.zeros(len(ynan))
lower = np.zeros(len(ynan))

for t in range(len(ynew)): 
    y_est_pred_sm_nan[t] = Z @ ks_res_nan.alpha_sm[:,:,t]
    upper[t] = y_est_pred_sm_nan[t] + ks_res_nan.eps_var[t,]
    lower[t] = y_est_pred_sm_nan[t] - ks_res_nan.eps_var[t,]
```


```python
plt.plot(u,ynan,label="Actual Data")
plt.plot(u,y_est_pred_sm_nan,c='yellow',label="Smotther Pred")
plt.plot(u,y_pred_nan,c='lawngreen',label="Filter Pred")
plt.axvline(x=u[n],c='r',label="Switch")
plt.legend()
plt.xlabel('Year')
plt.ylabel('GMSL deviation (mm)')
plt.show()
```

## Appendix A. Why didn't we learn the trend noise variance as well?
In the assignment above we have fixed $\sigma_\mu$ to a small value. Conceptually it would have been straightforward to learn also this parameter with the EM algorithm. However, unfortunately, the maximum likelihood estimate of $\sigma_\mu$ often ends up being too large to result in accurate _long term predictions_. The reason for this issue is that the structural model

\begin{align}
    y_t = \mu_t + \gamma_t + \varepsilon_t
\end{align}

is not a perfect description of reality. As a consequence, when learning the parameters the mismatch between the model and the data is compensated for by increasing the noise variances. This results in a trend component which does not only capture the long term trends of the data, but also seemingly random variations due to a model misspecification, possibly resulting in poor _long range predictions_.

Kitagawa (Introduction to Time Series Modeling, CRC Press, 2010, Section 12.3) discusses this issue and proposes two solutions. The first is a simple and pragmatic one: simply fix $\sigma_\mu^2$ to a value smaller than the maximum likelihood estimate. This is the approach we have taken in this assignment. The issue is of course that in practice it is hard to know what value to pick, which boild down to manual trial and error (or, if you are lucky, the designer of the lab assignment will tell you which value to use!).

The second, more principled, solution proposed by Kitagawa is to augment the model with a stationary AR component as well. That is, we model

\begin{align}
    y_t = \mu_t + \gamma_t + \nu_t + \varepsilon_t
\end{align}

where $\nu_t \sim$ AR$(p)$. By doing so, the stationary AR component can compensate for the discrepancies between the original structural model and the "true data generating process". It is straightforward to include this new component in the state space representation (how?) and to run the Kalman filter and smoother on the resulting model. Indeed, this is one of the beauties with working with the state space representation of time series data! However, the M-step of the EM algorithm becomes a bit more involved if we want to use the method to estimate also the AR coefficients of the $\nu$-component, which is beyond the scope of this lab assignment.

70ce47f2d632c865810959ef6e4aa3f208c247a1
ipynb
Jupyter Notebook
Lab2/tssl_lab2_Vysh.ipynb
vyspi394/Time-Series-Sequence-Learning
f9791b4bf36d07d2ba8b2f3c4649f851eded47ff
MIT
Lab2/tssl_lab2_Vysh.ipynb
vyspi394/Time-Series-Sequence-Learning
f9791b4bf36d07d2ba8b2f3c4649f851eded47ff
MIT
Lab2/tssl_lab2_Vysh.ipynb
vyspi394/Time-Series-Sequence-Learning
f9791b4bf36d07d2ba8b2f3c4649f851eded47ff
MIT
Qwen/Qwen-72B
1. YES
2. YES