# 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
from scipy import  linalg
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(u,y)
plt.xlabel('Year')
plt.ylabel('GMSL deviation (mm)')
plt.title("Data")
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:**

Dimensions is given as  Row * Column and it is given that d = dim($\alpha_{t}$) and $d_\eta = \dim(\eta_t)$

T = d * d

R = $d_\eta$ * 2

Z = 1 * d


**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
T_mu = np.array([[2,-1],[1,0]])
R_mu = np.array([[1],[0]])
Z_mu = np.array([[1 ,0 ]])
```


```python
print(T_mu.shape)
print(R_mu.shape)
print(Z_mu.shape)
```

    (2, 2)
    (2, 1)
    (1, 2)


**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
#s = 37 manual check
#coverting sample to year unit only
sampleYearly = u.astype(int)
unique_elements, counts_elements = np.unique(sampleYearly, return_counts=True)
#because of presence of last year count is decreasing by 1 so removing it as avg value is coming as 36 while without that it is 37 so going ahead with 37
#s = round(np.mean(np.delete(counts_elements , -1)))
s = round(np.mean(counts_elements))
# Removed assumptions 
print(f"Average number of observations per (whole) year and round to the closest integer is {s}")
# print(f"Year   : {unique_elements}")
# print(f"count per year : {counts_elements} ")
# print( f"avg of count {s}")
```

    Average number of observations per (whole) year and round to the closest integer is 36.0




**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:** For S as 36 seasonal component, state dimension will be S-1 * S-1 i.e 35*35


**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 = -1 * np.ones((1,35) , dtype=int)
T_gamma = np.append(T_gamma , np.append(np.identity(34) , np.zeros((34 ,1)) , axis = 1 ) , axis= 0)
R_gamma = np.zeros((35,1) , dtype= int)
R_gamma[0,0] = 1
Z_gamma = np.zeros((1,35) , dtype = int)
Z_gamma[0,0] = 1 
```


```python
print(T_gamma.shape)
print(R_gamma.shape)
print(Z_gamma.shape)
```

    (35, 35)
    (35, 1)
    (1, 35)


**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 = scipy.linalg.block_diag(T_mu , T_gamma)
R = scipy.linalg.block_diag(R_mu , R_gamma)
Z = np.concatenate((Z_mu , Z_gamma) , axis = 1 )
```


```python
print(f"Shape for T is {T.shape}")
print(f"Shape for R is {R.shape}")
print(f"Shape for Z is {Z.shape}")
```

    Shape for T is (37, 37)
    Shape for R is (37, 2)
    Shape for Z is (1, 37)



```python
#Z
```

Looking at above shape we can say that they corresponds to Q0 as dimensions are same as expected.

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
```

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 = y[0,]
# a1 = np.append(a1.reshape(1,1) , np.zeros((36,1)) , axis=0)
```


```python
a1 = np.zeros(37).reshape(37,1)
a1[0] = y[0] 
```


```python
P1 = 100 * np.identity(n = 37)
```


```python
P1.shape
```




    (37, 37)



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  # Module available in LISAM. kfs_res is a container class for storing the result.

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)) #renamed 
    y_pred = np.zeros(n)
    F_pred = np.zeros(n)
    
    #print(alpha_filt.shape)

    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
            #@ > dot product
            alpha_pred[:,:,t] = T @ alpha_filt[:,:,(t-1)] 
            #alpha_pred[:,:,t] = alpha_filt[:,:,(t-1)] 
            P_pred[:,:,t] = (T @ P_filt[:,:,(t-1)]) @ np.transpose(T) + R @ Q @ np.transpose(R)
        # Compute prediction of current output
        y_pred[t] = Z @ alpha_pred[:,:,t]
        F_pred[t] = ((Z @ P_pred[:,:,t]) @ np.transpose(Z)) +  H
        #F_pred[t] = Z @ P_pred[:,:,t] @ np.transpose(Z) 

        # 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_Gain = (P_pred[:,:,t] @ np.transpose(Z))/F_pred[t] 
          alpha_filt[:,:,t] = alpha_pred[:,:,t] + (K_Gain * (y[t] - y_pred[t]))
          #print(alpha_filt.shape)
          P_filt[:,:,t] = (np.identity(d) - (K_Gain @ Z)) @ P_pred[:,:,t]

    # Container for storing all results
    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
#it looks like we have some issue with n consideration as in Q7 it is asked to use 800 while in Q8 it is asked to use 700 
#Going with Q7 way with n as 800 datapoints 
n = 800
m = ndata-n
yTrain = y.copy()
#implemented such that for we can use it for long range prediction(i.e. adding Nan for future values) 
yTrain[801:] = np.nan
kf = kalman_filter(yTrain , model)

sd1 = kf.y_pred - np.sqrt(kf.F_pred)
sd2 = kf.y_pred + np.sqrt(kf.F_pred)
```


```python
#Data Plot 
plt.plot(u,y ,label = "Original Data")
plt.plot(u , kf.y_pred , label = "Kalman Pred")
#plt.scatter(u , sd1 , label = "Standard Deviation" , color = "green")
#plt.scatter(u , sd2 , label = "Standard Deviation" , color = "green")
plt.plot(u , sd1 , label = "Standard Deviation" , color = "green")
plt.plot(u , sd2 , label = "Standard Deviation" , color = "green")
plt.axvline(x=u[n])
plt.xlabel("Time")
plt.ylabel("GMSL")
plt.title("Kalman Filter")
plt.legend()
plt.show()
```

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


```python
#logLikLi = (-1/2)*(np.sum(np.log(np.absolute(kf.F_pred[:801,])) + (np.transpose(yTrain[:801,] - kf.y_pred[:801,]) / kf.F_pred[:801,]) @ (yTrain[:801,] - kf.y_pred[:801,])))
# -2085686.7746347892 old value kept for comparision
logLikLi = (-1/2)*(np.sum(np.log(np.abs(kf.F_pred[:801,])) + (yTrain[:801,] - kf.y_pred[:801,])**2 / kf.F_pred[:801,]))
print(f"Log-Likelihood for Training Data is {logLikLi}")

```

    Log-Likelihood for Training Data is -3199.819037286188


## 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 | \theta_{known}$)




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

_Hint: Look at Exercise Session 2_

**A:**

$\sigma_{\epsilon} ^ 2 $ = 1/n * $\sum_{t =1}^{n}$ ($\hat{\epsilon_{t}^2} + Var(\epsilon_{t} |Y_{n})  $)


$\sigma_{\eta} ^ 2$ = 1/n * $\sum_{t =1}^{n}$trace(($\hat{\eta_{t}^2}$ + Cov($\eta_t|Y_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
yTrain = y.copy()
yTrain = yTrain[:801]
# 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

def compute_theta(ks):
    """Implements the M-step of the EM-algorithm, based on the Kalman smoother results (ks). See Q11."""
    n = len(ks.eps_hat)
    sigma_epsilon2  = (1/n) * np.sum(ks.eps_hat**2 + ks.eps_var) #observation 
    #need to be calculated along diagonal only
    #sigma_neta2 = (1/(n-1)) * np.sum(ks.eta_hat**2 - ks.eta_cov) #seasonal 
    
    traceMatSum = 0
    #t_Mat = (ks.eta_hat.T @  ks.eta_hat + ks.eta_cov)    # dimensinal error because of grouping
    #  traceMatSum = np.trace(t_Mat)
    for t in range(n):
      t_Mat = (ks.eta_hat[:,:,t] @  ks.eta_hat[:,:,t].T + ks.eta_cov[:,:,t])  
      traceMatSum = traceMatSum + t_Mat[1,1]
   
    
    sigma_neta2 = (1/n) * traceMatSum
    theta_new = (sigma_epsilon2 , sigma_neta2)
    return theta_new


minChange = 0.001 #custom threshold set, could be changed based on expected diff 
logLikeLi = [] 
sigma_eps_Arr = []
sigma_seas_Arr = []
#Q = np.array([[sigma_trend**2, 0.], [0., sigma_seas**2]])  # Process noise covariance matrix #from above

sigma_eps2 = 1**2 #given

#changing r to 0 in-order to python indexing 
for r in range(0,num_iter):
    print(f"Number of iteration is : {r}")
    # E-step
    model = LGSS(T, R , Q, Z , sigma_eps2 , a1 , P1) #model Parameter
    kf = kalman_filter(yTrain , model)#Call Kalman filter with updated model 
    ks = kalman_smoother(yTrain, model, kf)
    
    #calculating log likelihood while ignoring constant as it is same accross all 
        
    
    logLikeLi.append((-1/2)*(np.sum(np.log(np.abs(kf.F_pred)) + (yTrain - kf.y_pred)**2 / kf.F_pred)))
    
    if r > 0:
      if (abs(logLikeLi[r] - logLikeLi[(r-1)]) > minChange):
        print(f"LogLikeLihood is {logLikeLi[r]}")
      else:
        print("Min change does not occur thus converging and loop break is applied")
        break;
    
    # M-step
    sigma_eps2 , sigma_seas2 = compute_theta(ks)
    sigma_eps_Arr.append(np.sqrt(sigma_eps2))
    sigma_seas_Arr.append(np.sqrt(sigma_seas2))
    Q = np.array([[sigma_trend**2, 0.], [0., sigma_seas2]])
    
```

    Number of iteration is : 0
    Number of iteration is : 1
    LogLikeLihood is -1866.8309059076453
    Number of iteration is : 2
    LogLikeLihood is -1737.1008966208346
    Number of iteration is : 3
    LogLikeLihood is -1705.5712890569646
    Number of iteration is : 4
    LogLikeLihood is -1693.6708986444391
    Number of iteration is : 5
    LogLikeLihood is -1686.6694873942026
    Number of iteration is : 6
    LogLikeLihood is -1681.1852264737176
    Number of iteration is : 7
    LogLikeLihood is -1676.4161295286615
    Number of iteration is : 8
    LogLikeLihood is -1672.1401988475586
    Number of iteration is : 9
    LogLikeLihood is -1668.2648409078197
    Number of iteration is : 10
    LogLikeLihood is -1664.7316664581126
    Number of iteration is : 11
    LogLikeLihood is -1661.4961324673172
    Number of iteration is : 12
    LogLikeLihood is -1658.5220925431386
    Number of iteration is : 13
    LogLikeLihood is -1655.7795117158917
    Number of iteration is : 14
    LogLikeLihood is -1653.2430877948818
    Number of iteration is : 15
    LogLikeLihood is -1650.8912800822036
    Number of iteration is : 16
    LogLikeLihood is -1648.7055856072982
    Number of iteration is : 17
    LogLikeLihood is -1646.669985593203
    Number of iteration is : 18
    LogLikeLihood is -1644.7705149168219
    Number of iteration is : 19
    LogLikeLihood is -1642.9949227213378
    Number of iteration is : 20
    LogLikeLihood is -1641.3324017506982
    Number of iteration is : 21
    LogLikeLihood is -1639.7733702059081
    Number of iteration is : 22
    LogLikeLihood is -1638.3092942185172
    Number of iteration is : 23
    LogLikeLihood is -1636.932542070822
    Number of iteration is : 24
    LogLikeLihood is -1635.6362634759316
    Number of iteration is : 25
    LogLikeLihood is -1634.4142888217982
    Number of iteration is : 26
    LogLikeLihood is -1633.2610444590969
    Number of iteration is : 27
    LogLikeLihood is -1632.1714809890707
    Number of iteration is : 28
    LogLikeLihood is -1631.141012167675
    Number of iteration is : 29
    LogLikeLihood is -1630.1654625446358
    Number of iteration is : 30
    LogLikeLihood is -1629.2410223401075
    Number of iteration is : 31
    LogLikeLihood is -1628.3642083589943
    Number of iteration is : 32
    LogLikeLihood is -1627.5318299741548
    Number of iteration is : 33
    LogLikeLihood is -1626.740959390622
    Number of iteration is : 34
    LogLikeLihood is -1625.988905546466
    Number of iteration is : 35
    LogLikeLihood is -1625.2731911189207
    Number of iteration is : 36
    LogLikeLihood is -1624.5915321959878
    Number of iteration is : 37
    LogLikeLihood is -1623.9418202463671
    Number of iteration is : 38
    LogLikeLihood is -1623.3221060799156
    Number of iteration is : 39
    LogLikeLihood is -1622.730585539519
    Number of iteration is : 40
    LogLikeLihood is -1622.1655867042268
    Number of iteration is : 41
    LogLikeLihood is -1621.6255584166215
    Number of iteration is : 42
    LogLikeLihood is -1621.1090599742745
    Number of iteration is : 43
    LogLikeLihood is -1620.6147518477114
    Number of iteration is : 44
    LogLikeLihood is -1620.141387306021
    Number of iteration is : 45
    LogLikeLihood is -1619.6878048476037
    Number of iteration is : 46
    LogLikeLihood is -1619.2529213464259
    Number of iteration is : 47
    LogLikeLihood is -1618.8357258361204
    Number of iteration is : 48
    LogLikeLihood is -1618.4352738634138
    Number of iteration is : 49
    LogLikeLihood is -1618.0506823515466
    Number of iteration is : 50
    LogLikeLihood is -1617.6811249204864
    Number of iteration is : 51
    LogLikeLihood is -1617.3258276180336
    Number of iteration is : 52
    LogLikeLihood is -1616.9840650203664
    Number of iteration is : 53
    LogLikeLihood is -1616.6551566657492
    Number of iteration is : 54
    LogLikeLihood is -1616.3384637891636
    Number of iteration is : 55
    LogLikeLihood is -1616.033386328767
    Number of iteration is : 56
    LogLikeLihood is -1615.7393601789618
    Number of iteration is : 57
    LogLikeLihood is -1615.4558546665437
    Number of iteration is : 58
    LogLikeLihood is -1615.1823702299414
    Number of iteration is : 59
    LogLikeLihood is -1614.9184362827227
    Number of iteration is : 60
    LogLikeLihood is -1614.66360924502
    Number of iteration is : 61
    LogLikeLihood is -1614.417470727909
    Number of iteration is : 62
    LogLikeLihood is -1614.179625857425
    Number of iteration is : 63
    LogLikeLihood is -1613.9497017259837
    Number of iteration is : 64
    LogLikeLihood is -1613.7273459603666
    Number of iteration is : 65
    LogLikeLihood is -1613.5122253964896
    Number of iteration is : 66
    LogLikeLihood is -1613.3040248517893
    Number of iteration is : 67
    LogLikeLihood is -1613.1024459873001
    Number of iteration is : 68
    LogLikeLihood is -1612.9072062520027
    Number of iteration is : 69
    LogLikeLihood is -1612.718037902721
    Number of iteration is : 70
    LogLikeLihood is -1612.534687093556
    Number of iteration is : 71
    LogLikeLihood is -1612.3569130292942
    Number of iteration is : 72
    LogLikeLihood is -1612.1844871776377
    Number of iteration is : 73
    LogLikeLihood is -1612.0171925358984
    Number of iteration is : 74
    LogLikeLihood is -1611.8548229477094
    Number of iteration is : 75
    LogLikeLihood is -1611.6971824659572
    Number of iteration is : 76
    LogLikeLihood is -1611.5440847585412
    Number of iteration is : 77
    LogLikeLihood is -1611.3953525536604
    Number of iteration is : 78
    LogLikeLihood is -1611.2508171216207
    Number of iteration is : 79
    LogLikeLihood is -1611.1103177905493
    Number of iteration is : 80
    LogLikeLihood is -1610.9737014934562
    Number of iteration is : 81
    LogLikeLihood is -1610.8408223444826
    Number of iteration is : 82
    LogLikeLihood is -1610.7115412420712
    Number of iteration is : 83
    LogLikeLihood is -1610.5857254972204
    Number of iteration is : 84
    LogLikeLihood is -1610.463248485114
    Number of iteration is : 85
    LogLikeLihood is -1610.3439893182465
    Number of iteration is : 86
    LogLikeLihood is -1610.2278325399295
    Number of iteration is : 87
    LogLikeLihood is -1610.1146678363139
    Number of iteration is : 88
    LogLikeLihood is -1610.0043897659975
    Number of iteration is : 89
    LogLikeLihood is -1609.896897505997
    Number of iteration is : 90
    LogLikeLihood is -1609.7920946126192
    Number of iteration is : 91
    LogLikeLihood is -1609.6898887968396
    Number of iteration is : 92
    LogLikeLihood is -1609.590191712637
    Number of iteration is : 93
    LogLikeLihood is -1609.4929187578161
    Number of iteration is : 94
    LogLikeLihood is -1609.3979888862357
    Number of iteration is : 95
    LogLikeLihood is -1609.3053244310327
    Number of iteration is : 96
    LogLikeLihood is -1609.2148509376677
    Number of iteration is : 97
    LogLikeLihood is -1609.1264970066754
    Number of iteration is : 98
    LogLikeLihood is -1609.0401941450857
    Number of iteration is : 99
    LogLikeLihood is -1608.9558766261744



```python
# for t in range(5):
#   print((ks.eta_hat[:,:,t] @  ks.eta_hat[:,:,t].T + ks.eta_cov[:,:,t]).shape)  
# #ks.eta_cov.shape
```


```python
plt.plot(range(0,len(sigma_eps_Arr)) , (sigma_eps_Arr) , label = "sigma_eps")
plt.plot(range(0,len(sigma_seas_Arr)) , (sigma_seas_Arr) , label = "sigma_eta")
plt.xlabel("Index")
plt.ylabel("Q")
plt.title("Traces of the parameter estimates")
plt.legend()
plt.show()
```


```python
#likelihood plot

plt.plot(range(0,len(logLikeLi)) , logLikeLi , label = "Loglikelihood")
plt.xlabel("Iteration")
plt.ylabel("Log Likelihood")
plt.title("Log Likelihood Vs Number of Iteration")
plt.legend()
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
#create model with new parameters 
#since last values will be used to direectly creating model with last Var values 
model = LGSS(T, R , Q, Z , sigma_eps2 , a1 , P1)
n = 800
m = ndata-n
yTrain = y.copy()
#implemented such that for we can use it for long range prediction(i.e. adding Nan for future values) 
yTrain[801:] = np.nan
kf = kalman_filter(yTrain , model)

sd1 = kf.y_pred - np.sqrt(kf.F_pred)
sd2 = kf.y_pred + np.sqrt(kf.F_pred)

```


```python
#Data Plot 
plt.plot(u,y ,label = "Original Data")
plt.plot(u , kf.y_pred , label = "Kalman Pred with EM Parameter")
#plt.scatter(u , sd1 , label = "Standard Deviation" , color = "green")
#plt.scatter(u , sd2 , label = "Standard Deviation" , color = "green")
plt.plot(u , sd1 , label = "Standard Deviation" , color = "green" , ls = "-.")
plt.plot(u , sd2 , label = "Standard Deviation" , color = "green" , ls = "-.")
plt.axvline(x=u[n])
plt.xlabel("Time")
plt.ylabel("GMSL")
plt.title("Kalman filter for long range prediction")
plt.legend()
plt.show()
```


```python
#Traning Data Log Likelihood
logLikLi = (-1/2)*(np.sum(np.log(np.abs(kf.F_pred[:801,])) + (yTrain[:801,] - kf.y_pred[:801,])**2 / kf.F_pred[:801,]))
print(f"Log-Likelihood for Training Data is {logLikLi}")
```

    Log-Likelihood for Training Data is -1608.8734813570054


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
ksSE = kalman_smoother(yTrain , model , kf) 
```


```python
#since s_t = z * alpha_sum 

n = len(yTrain)
s_tMean = []
s_tVar = []
for t in range(0,n):
  s_tMean.append(Z @ ksSE.alpha_sm[:,:,t])
  s_tVar.append(((Z @ ksSE.V[:,:,t]) @ np.transpose(Z)) + sigma_eps)

from pandas.core.common import flatten
s_tMean = list(flatten(s_tMean))
s_tVar = list(flatten(s_tVar))

sd1KS = s_tMean - np.sqrt(s_tVar)
sd2KS = s_tMean + np.sqrt(s_tVar)
```


```python
#Data Plot 
from pandas.core.common import flatten
s_tMean = list(flatten(s_tMean))
s_tVar = list(flatten(s_tVar))

sd1KS = s_tMean - np.sqrt(s_tVar)
sd2KS = s_tMean + np.sqrt(s_tVar)

plt.plot(u,y ,label = "Original Data")
plt.plot(u , s_tMean , label = "Kalman Smoothing Mean")
plt.scatter(u , sd1KS , label = "Standard Deviation" , color = "green" , ls = ":")
plt.scatter(u , sd2KS , label = "Standard Deviation" , color = "green" , ls = ":")
plt.axvline(x=u[801])
plt.xlabel("Time")
plt.ylabel("GMSL")
plt.title("Kalman smoother for computed smoothed estimates")
plt.legend()
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:**
To predict Kalman filter uses data till time t while in case of smoother we use all available (i.e. complete data).

**Training Data Points**

It can be seen that for Kalman filter prediction are changing rapidly since it only uses t-1 data for prediction while Kalman Smoother prediction is smoother since it is using future values to correct error in filtering prediction. 

**Validation data points**

Now while predicting validation data points, all future values are Nan so both predictions are similar as both access to the same data points only which is past values.


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
#ksSE.alpha_sm.shape
#ksSE.alpha_sm[0:2,0,:].sum(axis = 0).shape

plt.scatter(u , ksSE.alpha_sm[0:2,0,:].sum(axis = 0) , label = "Tread" )
plt.scatter(u , ksSE.alpha_sm[2:37,0,:].sum(axis = 0) , label = "Seasonal" )
plt.xlabel("Time")
plt.ylabel("GMSL")
plt.title("smoothed estimates signal components")
plt.legend()
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).


```python
yNew = y.copy()
yNew[300:400] = np.nan
```


```python
#Filter 
kfNew = kalman_filter(yNew , model)

#Smoother
ksNew = kalman_smoother(yNew , model , kfNew) 

n = len(yNew)
s_tMeanNew = []

for t in range(0,n):
  s_tMeanNew.append(Z @ ksSE.alpha_sm[:,:,t])

from pandas.core.common import flatten
s_tMeanNew = list(flatten(s_tMeanNew))
```


```python
plt.plot(u,y ,label = "Original Data")
plt.plot(u , kfNew.y_pred , label = "Kalman Filter")
plt.plot(u , s_tMeanNew , label = "Kalman Smoother")
plt.axvline(x=u[300])
plt.axvline(x=u[400])
plt.xlabel("Time")
plt.ylabel("GMSL")
plt.title("Kalman Filter and Kalman Smoother with missing data")
plt.legend()
plt.show()
```


```python
plt.plot(u[290:410],y[290:410] ,label = "Original Data")
plt.plot(u[290:410] , kfNew.y_pred[290:410] , label = "Kalman Filter")
plt.plot(u[290:410] , s_tMeanNew[290:410] , label = "Kalman Smoother")
plt.axvline(x=u[300] , color = "black")
plt.axvline(x=u[400] , color = "black")
plt.xlabel("Time")
plt.ylabel("GMSL")
plt.title("Kalman Filter and Kalman Smoother with missing data with subsection view")
plt.legend()
plt.show()
```

It has had been seen while predicting missing values, Filter predictions are rapidly changing which can predict changing amplitude, but after a point, it is not able to follow data properly.

Now smoother can follow trend and seasonality in data comparatively better when we have access better amount of past and future datapoints.

## 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.

62d0bee30e15f21a795ba0314a75a807bc224ff0
ipynb
Jupyter Notebook
Task2/tssl_lab2.ipynb
amannayak/Time-Series-and-Sequence-Analysis
543b0df70fe603c0ae08e9be8de632895f1bd24c
MIT
Task2/tssl_lab2.ipynb
amannayak/Time-Series-and-Sequence-Analysis
543b0df70fe603c0ae08e9be8de632895f1bd24c
MIT
Task2/tssl_lab2.ipynb
amannayak/Time-Series-and-Sequence-Analysis
543b0df70fe603c0ae08e9be8de632895f1bd24c
MIT
Qwen/Qwen-72B
1. YES
2. YES