# TSSL Lab 1 - Autoregressive models

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


```python
import pandas as pd # Loading data / handling data frames
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.pyplot import figure
from sklearn import linear_model as lm  # Used for solving linear regression problems
from sklearn.neural_network import MLPRegressor # Used for NAR model

from tssltools_lab1 import acf, acfplot # Module available in LISAM - Used for plotting ACF
```

## 1.1 Loading, plotting and detrending data

In this lab we will build autoregressive models for a data set corresponding to the Global Mean Sea Level (GMSL) over the past few decades. The data is taken from https://climate.nasa.gov/vital-signs/sea-level/ and is available on LISAM in the file `sealevel.csv`.

**Q1**: Load the data and plot the GMSL versus time. How many observations are there in total in this data set?

_Hint:_ With pandas you can use the function `pandas.read_csv` to read the csv file into a data frame. Plotting the time series can be done using `pyplot`. Note that the sea level data is stored in the 'GMSL' column and the time when each data point was recorded is stored in the column 'Year'.

**A1**:


```python
seaLevelData = pd.read_csv("sealevel.csv")
```


```python
#plot of GMSL vs Time 

plt.plot(seaLevelData['Year'] ,  seaLevelData['GMSL'])
plt.xlabel("Time")
plt.ylabel("GMSL")
plt.suptitle("GMSL VS Time")
plt.show()
```


```python
#total number of observation 

nObs = seaLevelData.shape[0]

print(f"Total number of obervation in given data is {nObs}")
```

    Total number of obervation in given data is 997


**Q2**: The data has a clear upward trend. Before fitting an AR model to this data need to remove this trend. Explain, using one or two sentences, why this is necessary.

**A2:** Since given data have a clear upward trend which introduces changing mean over time. Now to apply the AR model we need data with constant mean. Thus we have to detrend the data so stationary dataset can achieve and use for AR modeling and another benefit is that it enables discovering other aspects of the data like the autocorrelations and dependencies of the variables.

**Q3** Detrend the data following these steps:
1. Fit a straight line, $\mu_t=\theta_0 + \theta_1 u_t $ to the data based on the method of least squares. Here, $u_t$ is the time point when obervation $t$ was recorded.

    _Hint:_ You can use `lm.LinearRegression().fit(...)` from scikit-learn. Note that the inputs need to be passed as a 2D array.

    Before going on to the next step, plot your fitted line and the data in one figure.


2. Subtract the fitted line from $y_t$ for the whole data series and plot the deviations from the straight line.

**From now, we will use the detrended data in all parts of the lab.**

_Note:_ The GMSL data is recorded at regular time intervals, so that $u_{t+1} - u_t = $ const. Therefore, you can just as well use $t$ directly in the linear regression function if you prefer, $\mu_t=\theta_0 + \theta_1 t $.

**A3:**


```python
x = seaLevelData['Year'].to_numpy().reshape((-1 ,1))
y = seaLevelData['GMSL'].to_numpy()

model_lm = lm.LinearRegression().fit(x , y)

```


```python
theta0 = model_lm.intercept_
theta1 = model_lm.coef_
#Y_Pred = theta0 + theta1 * x
Y_Pred = model_lm.predict(x)

plt.plot(seaLevelData['Year'] ,  seaLevelData['GMSL'] , label = "Original Data")
plt.plot(x, Y_Pred, color='red' , label = "Fitted Line")
plt.xlabel("Time")
plt.ylabel("GMSL")
plt.legend()
plt.show()
```


```python
detrededData = y - Y_Pred

plt.plot(x, Y_Pred, color='red' , label = "Fitted Line")
plt.plot(x, detrededData , color='blue' , label = "Detrended Data")
plt.xlabel("Time")
plt.ylabel("GMSL")
plt.suptitle("Detrended Data GMSL Vs Time")
plt.legend()
plt.show()
```

**Q4:** Split the (detrended) time series into training and validation sets. Use the values from the beginning up to the 700th time point (i.e. $y_t$ for $t=1$ to $t=700$) as your training data, and the rest of the values as your validation data. Plot the two data sets.

_Note:_ In the above, we have allowed ourselves to use all the available data (train + validation) when detrending. An alternative would be to use only the training data also when detrending the model. The latter approach is more suitable if, either:
* we view the linear detrending as part of the model choice. Perhaps we wish to compare different polynomial trend models, and evaluate their performance on the validation data, or
* we wish to use the second chunk of observations to estimate the performance of the final model on unseen data (in that case it is often referred to as "test data" instead of "validation data"), in which case we should not use these observations when fitting the model, including the detrending step.

In this laboration we consider the linear detrending as a predetermined preprocessing step and therefore allow ourselves to use the validation data when computing the linear trend.

**A4:**


```python
trainDataY = detrededData[0:700 , ]
ValidDataY = detrededData[700: , ]

trainDataX = seaLevelData['Year'][0:700]
ValidDataX = seaLevelData['Year'][700:]
```


```python
plt.plot(trainDataX, trainDataY, color='blue')
plt.xlabel("Time")
plt.ylabel("GMSL")
plt.suptitle("Training Data")
plt.show()


plt.plot(ValidDataX, ValidDataY, color='blue')
plt.xlabel("Time")
plt.ylabel("GMSL")
plt.suptitle("Validation Data")
plt.show()

```

## 1.2 Fit an autoregressive model
We will now fit an AR$(p)$ model to the training data for a given value of the model order $p$.

**Q5**: Create a function that fits an AR$(p)$ model for an arbitrary value of p. Use this function to fit a model of order $p=10$ to the training data and write out (or plot) the coefficients.

_Hint:_ Since fitting an AR model is essentially just a standard linear regression we can make use of `lm.LinearRegression().fit(...)` similarly to above. You may use the template below and simply fill in the missing code.

**A5:**


```python
def fit_ar(y, p):
    """
    Fits an AR(p) model. The loss function is the sum of squared errors from t=p+1 to t=n.

    :param y: array (n,), training data points
    :param p: int, AR model order
    :return theta: array (p,), learnt AR coefficients
    """

    # Number of training data points
    n = y.shape[0]
    
    # Construct the regression matrix
    Phi = np.zeros((n-p,p), dtype= int)

    for j in range(p):
        Phi[:,j] = y[j:n+j-p]    
    
    # Drop the first p values from the target vector y
    yy = y[p:]  # yy = (y_{t+p+1}, ..., y_n)

  
    # Here we use fit_intercept=False since we do not want to include an intercept term in the AR model
    regr = lm.LinearRegression(fit_intercept=False)
    regr.fit(Phi,yy)    

    return regr.coef_
```


```python
#Calculating theta_hat from training data 
#set p = 10 as AR(10) model and detrained train data as training data
theta_hat = fit_ar( y = trainDataY , p = 10 )
plt.plot(range(len(theta_hat)), theta_hat, color='blue')
plt.ylabel("Theta")
plt.xlabel("Index")
plt.suptitle("Theta Coef")
plt.show()
```

**Q6:** Next, write a function that computes the one-step-ahead prediction of your fitted model. 'One-step-ahead' here means that in order to predict $y_t$ at $t=t_0$, we use the actual values of $y_t$ for $t<t_0$ from the data. Use your function to compute the predictions for both *training data* and *validation data*. Plot the predictions together with the data (you can plot both training and validation data in the same figure). Also plot the *residuals*.

_Hint:_ It is enought to call the predict function once, for both training and validation data at the same time.

**A6:**


```python
def predict_ar_1step(theta, y_target):
    """Predicts the value y_t for t = p+1, ..., n, for an AR(p) model, based on the data in y_target using
    one-step-ahead prediction.

    :param theta: array (p,), AR coefficients, theta=(a1,a2,...,ap).
    :param y_target: array (n,), the data points used to compute the predictions.
    :return y_pred: array (n-p,), the one-step predictions (\hat y_{p+1}, ...., \hat y_n) 
    """

    n = len(y_target)
    p = len(theta)
    
    # Number of steps in prediction
    m = n-p
    y_pred = np.zeros(m)

    for i in range(m):
      # if i == 0 :
      #   prevY = y_target[i:p+i]
      # #print(i)
      # #print(prevY[i:p+i])
      # y_pred[i] = np.flip(prevY[i:p+i]).dot(theta)  
      # #print(y_pred[i])
      # prevY = np.append(prevY , y_pred[i])
      # #print(prevY[i:p+i])

      y_pred[i] = y_target[i:p+i].dot(theta)
    return y_pred
```


```python
YPred_Train = predict_ar_1step(theta = theta_hat , y_target = trainDataY)
YPred_Valid = predict_ar_1step(theta = theta_hat , y_target = ValidDataY)

```


```python
#plot for predicted values for train and validation data 
plt.plot(trainDataX , trainDataY , label = "Training Data")
plt.plot(trainDataX[len(theta_hat):,] , YPred_Train , label = "Prediction for Training Data")
plt.plot(ValidDataX , ValidDataY , label = "Validation Data")
plt.plot(ValidDataX[len(theta_hat):,] , YPred_Valid , label = "Prediction for Validation Data")
#figure(num=None, figsize=(8, 6), dpi=80, facecolor='w', edgecolor='k')
plt.xlabel("Time")
plt.ylabel("GMSL")
plt.legend()
plt.show()

```

**Q7:** Compute and plot the autocorrelation function (ACF) of the *residuals* only for the *validation data*. What conclusions can you draw from the ACF plot?

_Hint:_ You can use the function `acfplot` from the `tssltools` module, available on the course web page.

**A7:**


```python
help(acfplot)
```

    Help on function acfplot in module tssltools_lab1:
    
    acfplot(x, lags=None, conf=0.95)
        Plots the empirical autocorralation function.
        
        :param x: array (n,), sequence of data points
        :param lags: int, maximum lag to compute the ACF for. If None, this is set to n-1. Default is None.
        :param conf: float, number in the interval [0,1] which specifies the confidence level (based on a central limit
                     theorem under a white noise assumption) for two dashed lines drawn in the plot. Default is 0.95.
        :return:
    


The plot indicates that the sign of the auto-correlation changes depending on the value of the coefficients. Moreover, we can see that the decay in the ACF is relatively slow, meaning the process takes a long to forget its past.


```python
resiValidation = ValidDataY[10:] - YPred_Valid

AutoCoR = acfplot(x = resiValidation)


```

## 1.3 Model validation and order selection
Above we set the model order $p=10$ quite arbitrarily. In this section we will try to find an appropriate order by validation.

**Q8**: Write a loop in which AR-models of orders from $p=2$ to $p=150$ are fitted to the data above. Plot the training and validation mean-squared errors for the one-step-ahead predictions versus the model order.

Based on your results:
- What is the main difference between the changes in training error and validation error as the order increases? 
- Based on these results, which model order would you suggest to use and why?

_Note:_ There is no obvious "correct answer" to the second question, but you still need to pick an order an motivate your choice!



```python
MSE_Train = []
MSE_Valid = []

# def ModelCompare(y,p):
#   theta_hat = fit_ar(y, p)
#   predY = predict_ar_1step(theta = theta_hat , y_target = y)
#   residuals  = y[len(theta_hat):] - predY
#   mse = sum(residuals**2) / len(residuals)
#   return mse

for i in range(2,151):
  #calculate theta_hat
  theta_hat = fit_ar(trainDataY, i)
  #pred for train Data
  predYtrain = predict_ar_1step(theta = theta_hat , y_target = trainDataY)
  residualsYtrain  = trainDataY[len(theta_hat):] - predYtrain
  MSE_Train.append(sum(residualsYtrain**2) / len(residualsYtrain))

  #pred for Valid Data
  predYValid = predict_ar_1step(theta = theta_hat , y_target = ValidDataY)
  residualsYvalid  = ValidDataY[len(theta_hat):] - predYValid
  MSE_Valid.append(sum(residualsYvalid**2) / len(residualsYvalid))

```


```python
plt.plot(range(2,151) , MSE_Train , label = "MSE Train Data")
plt.plot(range(2,151) , MSE_Valid , label = "MSE Valid Data")
plt.xlabel("p")
plt.ylabel("MSE")
plt.suptitle("MSE vs Model Order(P)")
plt.legend()
plt.show()
```

**A8:**

Looking at the above plot of MSE vs Model Order P we can see that MSE for training data decreases with increased model order p while for validation data it is marginally changing between 4.6 and 5.

While around model order 107 both values decrease and post that MSE for Validation Data starts to increase.

Based on the above result we can consider model order as 10.



**Q9:** Based on the chosen model order, compute the residuals of the one-step-ahead predictions on the *validation data*. Plot the autocorrelation function of the residuals. What conclusions can you draw? Compare to the ACF plot generated above for p=10.


```python
#at model order of 110

#calculate theta_hat
theta_hat = fit_ar(trainDataY, 107)
#pred for Valid Data
predYValid = predict_ar_1step(theta = theta_hat , y_target = ValidDataY)
residualsValid107  = ValidDataY[len(theta_hat):] - predYValid

AutoCoR110 = acfplot(x = residualsValid107)

```

The plot suggests a decrease in the ACF compared to when p=10. Because we are using a much higher order now which means we
consider 107 previous times when predicting a new value. This leads to a
lower correlation as the lags used are far apart from the predictions
point.

## 1.4 Long-range predictions
So far we have only considered one-step-ahead predictions. However, in many practical applications it is of interest to use the model to predict further into the future. For intance, for the sea level data studied in this laboration, it is more interesting to predict the level one year from now, and not just 10 days ahead (10 days = 1 time step in this data).

**Q10**: 
Write a function that simulates the value of an AR($p$) model $m$ steps into the future, conditionally on an initial sequence of data points. Specifically, given $y_{1:n}$ with $n\geq p$ the function/code should predict the values

\begin{align}
    \hat y_{t|n} &= \mathbb{E}[y_{t} | y_{1:n}], & t&=n+1,\dots,n+m.
\end{align}

Use this to predict the values for the validation data ($y_{701:997}$) conditionally on the training data ($y_{1:700}$) and plot the result.

_Hint:_ Use the pseudo-code derived at the first pen-and-paper session.

**A10:**


```python
def simulate_ar(y, theta, m):
    """Simulates an AR(p) model for m steps, with initial condition given by the last p values of y
    
    :param y: array (n,) with n>=p. The last p values are used to initialize the simulation.
    :param theta: array (p,). AR model parameters,
    :param m: int, number of time steps to simulate the model for.
    """

    p = len(theta)    
    y_sim = np.zeros(m)
    phi = np.flip(y[-p:].copy()) # (y_{n-1}, ..., y_{n-p})^T - note that y[ntrain-1] is the last training data point
    
    #Last Element Index which is to be removed 
    #popOutIndex = len(phi) - 1  

    theta = np.flip(theta.copy())

    for i in range(m):
        y_sim[i] = phi.dot(theta)
        #phi = np.append(y_sim[i] , np.delete(phi , popOutIndex))      
        phi = np.append(y_sim[i] , np.delete(phi , -1))      
    return y_sim


```


```python
theta_Gen = fit_ar(y = trainDataY , p = 27)
m = 997 - 700

ySimPredValid = simulate_ar( y = trainDataY , theta = theta_Gen , m = m)

plt.plot(ValidDataX , ValidDataY , label = "Original Validation Data")
plt.plot(ValidDataX , ySimPredValid , label = "Simulated Validation Data")

#figure(num=None, figsize=(8, 6), dpi=80, facecolor='w', edgecolor='k')
plt.xlabel("Time")
plt.ylabel("GMSL")
plt.legend()
plt.show()

```

**Q11:** Using the same function as above, try to simulate the process for a large number of time steps (say, $m=2000$). You should see that the predicted values eventually converge to a constant prediction of zero. Is this something that you would expect to see in general? Explain the result.

**A11**






```python
ySimLargeM = simulate_ar( y = trainDataY , theta = theta_Gen , m = 2000)

plt.plot(range(2000) , ySimLargeM , label = "Simulated Data with m = 2000")
plt.xlabel("Time")
plt.ylabel("GMSL")
plt.legend()
plt.show()
```


```python
AutoCoR_YSimLarge = acfplot(x = ySimLargeM)
```

It can be seen that for such large prediction, predicted values end up being a function of Noise as after a point all past values considered are newly predicted values. Now Noise is distributed with mean = 0 thus it is expected to see a convergence of mean to 0.

It can be seen from the ACF plot that after a large lag like 1000, autocorrelation between values is very low and now they are only random walk generated by the distribution of noise function which has mean as 0, thus new predicted expected values are converging to 0. 

## 1.5 Nonlinear AR model
 In this part, we switch to a nonlinear autoregressive (NAR) model, which is based on a feedforward neural network. This means that in this model the recursive equation for making predictions is still in the form $\hat y_t=f_\theta(y_{t-1},...,y_{t-p})$, but this time $f$ is a nonlinear function learned by the neural network. Fortunately almost all of the work for implementing the neural network and training it is handled by the `scikit-learn` package with a few lines of code, and we just need to choose the right structure, and prepare the input-output data.   

**Q12**: Construct a NAR($p$) model with a feedforward (MLP) network, by using the `MLPRegressor` class from `scikit-learn`. Set $p$ to the same value as you chose for the linear AR model above. Initially, you can use an MLP with a single hidden layer consisting of 10 hidden neurons. 
Train it using the same training data as above and plot the one-step-ahead predictions as well as the residuals, on both the training and validation data. 

_Hint:_ You will need the methods `fit` and `predict` of `MLPRegressor`. Read the user guide of `scikit-learn` for more details. Recall that a NAR model is conceptuall very similar to an AR model, so you can reuse part of the code from above.

**A12:**


```python
def phiCal(y , p):
  n = y.shape[0]
  
  Phi = np.zeros((n-p,p), dtype= int)

  for j in range(p):
      Phi[:,j] = y[j:(n+j-p)]    
  return Phi
p = 108
yModel = trainDataY[p:]  

model = MLPRegressor(hidden_layer_sizes=(10,) , activation= "relu")
reg = model.fit(phiCal(y = trainDataY , p = 108), yModel)
n = trainDataY.shape[0]
dataValid_Y = np.append(trainDataY[(n - p): ] , ValidDataY)
predictedVal = reg.predict(phiCal(y = dataValid_Y , p = 108))
```

    /usr/local/lib/python3.6/dist-packages/sklearn/neural_network/_multilayer_perceptron.py:571: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (200) reached and the optimization hasn't converged yet.
      % self.max_iter, ConvergenceWarning)



```python
reg.coefs_[0].shape
```




    (108, 10)




```python
plt.plot(ValidDataX , predictedVal , label = "NAR Prediction with p = 108")
plt.plot(ValidDataX , ValidDataY , label = "Original Data")
plt.xlabel("Time")
plt.ylabel("GMSL")
plt.legend()
plt.show()
```


```python
NARresidual = ValidDataY - predictedVal
plt.plot(ValidDataX , ValidDataY , label = "Original Data")
plt.plot(ValidDataX , NARresidual , label = "Residuals")
plt.xlabel("Time")
plt.ylabel("GMSL")
plt.legend()
plt.show()
```

**Q13:** Try to expirement with different choices for the hyperparameters of the network (e.g. number of hidden layers and units per layer, activation function, etc.) and the optimizer (e.g. `solver` and `max_iter`).

Are you satisfied with the results? Why/why not? Discuss what the limitations of this approach might be.

**A13:**

We have tried with multiple settings, and with a higher number of hidden layers and other different parameters, we could see that residual for training data was significantly low while there was high variance for Validation data, which in turn suggested overfitting of the model. 

The below model captures non-linearity in data with a comparatively lower residual, thus performs better comparatively. 

If comparing the above two models below is satisfactory while it still has a high residual for validation data.


```python
model = MLPRegressor(hidden_layer_sizes=(30,10) , activation= "tanh" ,  
                     solver="sgd" , learning_rate = "adaptive" , max_iter = 200)
reg = model.fit(phiCal(y = trainDataY , p = 108), yModel)
n = trainDataY.shape[0]
#dataValid_Y = np.append(trainDataY[(n - p): ] , ValidDataY)
predictTrainData = reg.predict(phiCal(y = trainDataY , p = 108))
predictTrainDataResidual = trainDataY[p:] - predictTrainData
predictedVal = reg.predict(phiCal(y = dataValid_Y , p = 108))
NARresidualOpt = ValidDataY - predictedVal
plt.plot(ValidDataX , NARresidual , label = "Pre Changes Residuals A12 ")
plt.plot(ValidDataX , NARresidualOpt , label = "Post Changes Residuals A13")
plt.plot(trainDataX[p:] , predictTrainDataResidual , label = "Train Pred Residuals")
plt.xlabel("Time")
plt.ylabel("GMSL")
plt.legend()
plt.show()
```


```python
plt.plot(ValidDataX , NARresidual , label = "Pre Changes Residuals :: A12 ")
plt.plot(ValidDataX , NARresidualOpt , label = "Post Changes Residuals :: A13")
plt.xlabel("Time")
plt.ylabel("GMSL")
plt.legend()
plt.show()
```

30cd0e571beb2a053efc6f4d9e08658b394141b7
ipynb
Jupyter Notebook
Task1/tssl_lab1.ipynb
amannayak/Time-Series-and-Sequence-Analysis
543b0df70fe603c0ae08e9be8de632895f1bd24c
MIT
Task1/tssl_lab1.ipynb
amannayak/Time-Series-and-Sequence-Analysis
543b0df70fe603c0ae08e9be8de632895f1bd24c
MIT
Task1/tssl_lab1.ipynb
amannayak/Time-Series-and-Sequence-Analysis
543b0df70fe603c0ae08e9be8de632895f1bd24c
MIT
Qwen/Qwen-72B
1. YES
2. YES