# TSSL Lab 1 - Autoregressive models

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 matplotlib.pyplot as plt
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
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

gmsl = pd.read_csv('sealevel.csv')
print(f"Number of observations: {len(gmsl.index)}")

sns.lineplot(x='Year', y='GMSL', data=gmsl)
```

**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:**

Detrending the data means removing the global trend of the data over time (removing the mean function in this case). This is done by removing the mean function from the residuals. Benefit on this is that it enables discovering other aspects of the data like the auto correlations and dependences 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 = gmsl['Year'].to_numpy()
y = gmsl['GMSL'].to_numpy()

X = X.reshape(-1, 1)
#y = y.reshape(-1, 1)

model = lm.LinearRegression().fit(X,y)
```


```python
ypred = model.predict(X)

plt.plot(X, y,  color='blue')
plt.plot(X, ypred, color='red', linewidth=3)
```


```python
detrended = y - ypred

plt.plot(X,detrended, color='blue')
```

**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
# splitting the feature
Xtrain = X[0:700,]
Xvalid = X[700:,]

# splitting the target
detrendedtrain = detrended[0:700,]
detrendedvalid = detrended[700:,]
```

## 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 = len(y) # <COMPLETE THIS LINE>
    
   # Construct the regression matrix
    Phi = np.zeros((n-p,p), dtype= int)
 
     #print("phi")
     #print(Phi)
     #print(Phi.shape)
    for j in range(p):
        Phi[:,j] = y[p-j:n-j] # <COMPLETE THIS LINE>
    
    # 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
theta = fit_ar(detrendedtrain, 10)
theta
```




    array([ 1.04858020e+00,  1.03995740e-02,  6.87624028e-03,  6.22803308e-03,
            5.70699995e-03,  6.59110282e-03, -2.70953833e-03, -9.16897128e-03,
            6.66092353e-04, -7.84066761e-03])



**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)
    prevY = y_target[:10,]
    for i in range(m):
        # <COMPLETE THIS CODE BLOCK>
        #y_pred[i] = y_target[i+p:i].dot(theta) # <COMPLETE THIS LINE>
        prevY = prevY[i:p+i]
        y_pred[i] = prevY.dot(theta)  
        prevY = np.append(prevY , y_pred)

    return y_pred
```


```python
pred_train = predict_ar_1step(theta, detrendedtrain)
pred_valid = predict_ar_1step(theta, detrendedvalid)
```


```python
pred_valid.shape
```




    (287,)




```python
plt.plot(Xtrain[10:,], detrendedtrain[10:,],  color='blue')
plt.plot(Xtrain[10:,], pred_train, color='red')
```


```python
plt.plot(Xvalid[10:,], detrendedvalid[10:,],  color='blue')
plt.plot(Xvalid[10:,], pred_valid, color='red')
```


```python
pred_valid
```




    array([ -3.1360453 ,  -5.13169701,  -6.02434681,  -8.64803883,
            -3.32547717,  -6.40300014,  -5.43193911,  -9.14305106,
            -6.77966579,  -9.64969393,  -3.32853436,  -5.48956174,
            -6.41328221,  -9.17417627,  -3.62319396,  -6.85667756,
            -5.79996074,  -9.7194242 ,  -7.21849221, -10.18513903,
            -3.53454887,  -5.87249929,  -6.82839194,  -9.73422364,
            -3.94406179,  -7.34034545,  -6.19206856, -10.33103593,
            -7.68422883, -10.75111814,  -3.75511776,  -6.28227117,
            -7.27146231, -10.33038623,  -4.28970712,  -7.85589204,
            -6.60979393, -10.98001137,  -8.17852996, -11.3494737 ,
            -3.99134905,  -6.72076084,  -7.74439965, -10.96501006,
            -4.66186302,  -8.40532135,  -7.05476508, -11.6686053 ,
            -8.70315551, -11.98216617,  -4.24443569,  -7.18998185,
            -8.24923807, -11.64059071,  -5.0623759 ,  -8.99076027,
            -7.52871332, -12.39921028,  -9.25997793, -12.65128186,
            -4.51566188,  -7.69208663,  -8.78814786, -12.35978256,
            -5.49321264,  -9.61446611,  -8.03347973, -13.17436539,
            -9.85098976, -13.35904127,  -4.80640979,  -8.22937575,
            -9.36344427, -13.12540872,  -5.95646819, -10.27883465,
            -8.57102227, -13.99676538, -10.47831161, -14.10780804,
            -5.11816678,  -8.8043078 ,  -9.97759696, -13.94047146,
            -6.45437367, -10.98640864,  -9.14342339, -14.8692706 ,
           -11.14420071, -14.90009846,  -5.45253311,  -9.41951   ,
           -10.63323998, -14.8081634 ,  -6.98930491, -11.73988689,
            -9.75289808, -15.79491744, -11.85106004, -15.73859161,
            -5.81123021, -10.07778934, -11.3331824 , -15.73187928,
            -7.56379158, -12.54213388, -10.40180253, -16.77692951,
           -12.60144809, -16.62614019,  -6.19610949, -10.78214454,
           -12.08041963, -16.71522847,  -8.18052694, -13.39619005,
           -11.09264338, -17.81872952, -13.39808925, -17.56578205,
            -6.60916179, -11.53577882, -12.8781454 , -17.76204835,
            -8.84237807, -14.30528266, -11.82808751, -18.92395199,
           -14.24388494, -18.56075252,  -7.05252751, -12.34211334,
           -13.72976461, -18.87641837,  -9.55239687, -15.27283745,
           -12.61097257, -20.09645668, -15.14192546, -19.61449753,
            -7.52850734, -13.20480168, -14.63890683, -20.06267508,
           -10.31383173, -16.30249094, -13.44431819, -21.34034304,
           -16.09550266, -20.73068762,  -8.03957391, -14.12774512,
           -15.6094408 , -21.32542804, -11.13013988, -17.39810357,
           -14.33133793, -22.65996547, -17.10812351, -21.91323282,
            -8.58838404, -15.11510894, -16.64548974, -22.66957668,
           -12.00500062, -18.56377366, -15.27545203, -24.05994968,
           -18.18352449, -23.16629862,  -9.17779195, -16.17133983,
           -17.75144768, -24.10032831, -12.9423293 , -19.80385231,
           -16.2803011 , -25.54521004, -19.32568701, -24.49432291,
            -9.8108634 , -17.30118427, -18.93199682, -25.62321711,
           -13.94629223, -21.12295922, -17.34976062, -27.12096814,
           -20.5388539 , -25.90203403, -10.49089067, -18.50970825,
           -20.19212592, -27.24412436, -15.02132256, -22.52599949,
           -18.48795652, -28.79277249, -21.82754698, -27.39447011,
           -11.22140872, -19.80231811, -21.53714997, -28.96929993,
           -16.17213715, -24.0181816 , -19.69928172, -30.56651957,
           -23.19658577, -28.97699954, -12.00621232, -21.18478278,
           -22.97273094, -30.80538505, -17.40375451, -25.60503649,
           -20.98841387, -32.44847621, -24.65110753, -30.65534293,
           -12.84937438, -22.66325748, -24.50489999, -32.75943655,
           -18.72151388, -27.29243788, -22.36033428, -34.44530348,
           -26.19658864, -32.4355964 , -13.75526563, -24.2443088 ,
           -26.14008096, -34.83895252, -20.13109561, -29.0866239 ,
           -23.82034807, -36.56408212, -27.83886736, -34.32425648,
           -14.72857542, -25.9349415 , -27.88511544, -37.05189965,
           -21.63854273, -30.99422015, -25.37410576, -38.81233961,
           -29.58416821, -36.32824657, -15.77433414, -27.74262694,
           -29.7472893 , -39.40674219, -23.25028399, -33.02226426,
           -27.02762628, -41.19807908, -31.43912787, -38.4549452 ,
           -16.89793702, -29.6753334 , -31.73436104, -41.91247277,
           -24.97315838, -35.17823201, -28.78732155, -43.72981005,
           -33.41082294, -40.71221611, -18.10516957, -31.74155819,
           -33.85459179, -44.57864506, -26.81444116, -37.47006518,
           -30.66002273, -46.41658122, -35.50679944, -43.10844028,
           -19.40223475, -33.95036196, -36.11677726, -47.41540854,
           -28.78187166, -39.90620122, -32.65300827])




```python
plt.plot(Xvalid[700:,], detrendedvalid[700:,],  color='blue')
plt.plot(Xvalid[700:,], pred_valid, color='red', linewidth=3)
```

**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. 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:
    




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


**A8:**


```python

```



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

```



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

    for i in range(m):
        y_sim[i] = ... # <COMPLETE THIS LINE>
        # <COMPLETE THIS CODE BLOCK>
    
    return y_sim
    
```


```python

```

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

```

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

```


```python

```

**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:**


```python

```


```python

```

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