Local Linear Trend Model for time series analysis on Stan

Overview

On this article, I’ll make the local linear trend model to artificial time series data by Stan. Before, I made a local level model on Stan on the article, Local level model to time series data on Stan.
By adding the slope variable to that model, I’ll make the local linear trend model.
I used this book as reference.

Local Linear Trend Model

Actually, I’m not econometrics-oriented person and still on the way of studying how to analyze time series data. So, about local linear trend model, if something is wrong, please let me know.

Before, I made local level model and the model can be expressed by the following equation.
$u_t \sim Normal(u_{t-1}, σ_ξ^2)$
$y_t \sim Normal(u_t, σ_e^2$

$σ_e^2$ is observation disturbance. $σ_ξ^2$ is level disturbance.
On the local linear trend model, we need to add one more variable, slope variable.
$v_t \sim Normal(v_{t-1}, σ_ζ^2)$
$u_t \sim Normal(u_{t-1}+v_{t-1}, σ_ξ^2)$
$y_t \sim Normal(u_t, σ_e^2$

Data

The data is time series artificial one. Here, I set priority on using same data as other articles. So actually, here, don’t care about if the model is appropriate to the data or not.

import numpy as np

# make data
np.random.seed(59)
data = []
start_point = 10.0
data.append(start_point)
for i in range(100):
    if i == 0:
        temp = start_point
    else:
        next_temp = temp + np.random.normal(0, 1)
        data.append(next_temp)
        temp = next_temp

By plot, we can check how the data is.

import matplotlib.pyplot as plt
plt.plot(list(range(len(data))), data)
plt.show()

Stan Modeling

By the same manner as the local level model, I’ll write the model on Stan. When we compare with the local level model’s code, only difference is the existence of vector v.

data {
    int N;
    vector[N] X;
}

parameters {
    vector[N] u;
    vector[N] v;
    real<lower=0> s_u;
    real<lower=0> s_v;
    real<lower=0> s_x;
}

model {
    v[2:N] ~ normal(v[1:N-1], s_v);
    u[2:N] ~ normal(u[1:N-1] + v[1:N-1], s_u);
    X ~ normal(u, s_x);
}

After saving the code above to the file ‘local_linear_trend_model.stan’, we can execute on python.

Execute

On python, we can read the stan code and execute it.

import pystan

data_feed = {'X': data, 'N': len(data)}
fit = pystan.stan(file='local_linear_trend_model.stan', data=data_feed, iter=1000)

We can easily plot the outcome.

fit.plot()

To check the outcome, we can extract the sampled points. Here, I’ll check only $u$.

samples = fit.extract(permuted=True)
u_mean = samples['u'].mean(axis=0)

By plotting the data and the sampled $u$ at once, we can check well.

plt.plot(list(range(len(data))), data)
plt.plot(list(range(len(data))), u_mean)
plt.show()