# How to fit an autoregressive (AR(1)) model with trend and/or seasonality to a time series?

I want to test a model I have on a time series. The model is that the time series adapts to a trend $$f(t)$$ with a speed $$\alpha$$. There is also noise in the model. So, the time series is a function plus an AR(1) process. Mathematically the model is as follows:

$$y_{t + 1} = (1 - \alpha)y_{t} + \alpha f(t) + \epsilon_t$$ where $$\epsilon_t \sim \mathcal{N}(0,1)$$ and $$\alpha$$ is the parameter I want to fit. The only thing we know about $$f(t)$$ is that it changes slowly. I want to obtain $$\alpha$$, the speed at which the system adjusts to the external trend ($$f(t)$$).

Is there a way to do such fit?

My naive attempt was to use python library statsmodel to decompose the time series. After decomposing, I tried to fit an AR(1) model to the residuals to estimate $$\alpha$$ but I am not sure if this is mathematically founded. Below I show a toy example of my attempt using a very simple seasonality function $$f(t) = A(1 + \sin(t))$$, where $$A$$ is some amplitude.

The simulated time series.

Decomposition of time series. Residual is not the AR(1) process I hoped for.

Below the code, I used to generate the examples and the actual fit (which fails).

import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.tsa.arima_model import ARMA

#defining the trend function
def trend(t, amp=1):
return amp*(1 + np.sin(t))

#length of time series
n_time_steps = 250
#amplitud of time series
amplitud=10
#initializing the time series
time_series = np.zeros(n_time_steps)
time_series[0] = trend(0, amplitud)

#The AR(1) parameter. Our goal will be to find this parameter.
alpha = 0.1
#making the time series
for t in range(1,n_time_steps):
time_series[t] = (1 - alpha)*time_series[t - 1] + alpha*trend(t, amp=amplitud) + alpha*np.random.normal(0,1)

#passing the time series to a pandas format
dates = sm.tsa.datetools.dates_from_range('2000m1', length=len(time_series))
time_series_pd= pd.Series(time_series, index=dates)
#decomposing the time series
res = sm.tsa.seasonal_decompose(time_series_pd)
#fitting the AR(1) model
mod = ARMA(list(res.resid[6:-6]), order=(1,0))
ar1_fit = mod.fit()
#The alpha parameter is...
print(1 - ar1_fit.params[1]) #... wrong


I made the time series with a parameter $$\alpha = 0.1$$ but the fit was of $$0.4$$. I ran the code several times with different values of $$\alpha$$ and the resulting fit was always around 0.4 so the naive method is definitely not working.

My question is, is there a way in which we can fit $$\alpha$$ to an observed time series, where $$f(t)$$ can be a trend, seasonality or both?

I have very little knowledge of time series so references to similar problems would also be greatly appreciated.

Edits To answer the request in comments I add values below. Random seed is set to 1.

[10. , 10.72595885, 11.52358917, 11.34826863, 10.45888924, 9.48023122, 9.38369266, 10.07843437, 10.9787728 , 11.27459181, 10.63483361, 9.49017098, 8.87688341, 9.435671 , 10.57884086, 11.06298176, 10.57606807, 9.39747464, 8.72211285, 9.04041155, 10.22283847, 11.11969964, 11.13273664, 10.2063528 , 9.23146764, 9.13337707, 10.05086883, 10.93238285, 11.15283698, 10.28641578, 9.26554459, 8.73949121, 9.3985384 , 10.40135548, 10.75713964, 10.26221266, 9.39586217, 8.88760798, 9.2262743 , 10.26346529, 10.93825431, 10.67017454, 9.61173689, 8.72868546, 8.86633691, 9.64753428, 10.70219519, 10.63893072, 9.71554529, 8.87046796, 8.82872384, 9.67172944, 10.75581672, 10.91077525, 10.25928348, 9.19834758, 8.61558807, 9.25163173, 10.48970349, 11.10617916, 10.46977925, 9.5732737 , 8.85475124, 8.97903685, 10.00109005, 10.84209007, 10.75244881, 9.90575842, 9.04823385, 8.99914474, 9.84467462, 10.71211039, 10.80225592, 10.21593102, 9.19557713, 9.03864446, 9.83039095, 10.89768424, 11.41492796, 10.88613347, 9.74998291, 9.12889748, 9.77000435, 10.7379209 , 11.29238873, 11.01397586, 9.92263867, 9.0497853 , 9.06670674, 10.06797882, 10.91946614, 10.75447231, 9.8593913 , 8.77796022, 8.65757981, 9.35618397, 10.42261619, 10.81248526, 10.17407762, 9.16612511, 8.68497338, 9.25949543, 10.211867 , 10.82158327, 10.437084 , 9.33053699, 8.74546459, 9.04645056, 9.9856934 , 10.99536028, 10.8958873 , 9.94829915, 9.08498071, 9.0035942 , 9.91447592, 11.07908555, 11.20526544, 10.37895494, 9.37160344, 8.80010115, 9.46337901, 10.58708338, 10.78607294, 10.32424094, 9.16009977, 8.62251457, 8.87281301, 10.05402341, 10.72336744, 10.46214931, 9.51550051, 8.76102929, 8.72676343, 9.75099838, 10.81271358, 10.83078971, 9.92121378, 9.04353229, 8.84804249, 9.93844871, 10.95795999, 11.32510231, 10.58624608, 9.45974764, 8.87414693, 9.58640768, 10.598634 , 11.11499991, 10.76228583, 9.75698432, 8.98065554, 9.44689388, 10.42816038, 11.17091396, 11.00807861, 9.94053838, 9.1140563 , 8.9552398 , 9.88206657, 10.63387079, 10.92697227, 10.00724925, 8.91875369, 8.59912546, 9.34680328, 10.31753613, 10.64674645, 9.8817167 , 8.80257229, 8.23773322, 8.79169434, 9.89682222, 10.48888907, 10.25241797, 9.39400307, 8.7057619 , 8.99374914, 9.89105015, 10.74620886, 10.77312919, 9.91989239, 8.95413862, 8.87388559, 9.47038948, 10.45047886, 10.68004816, 9.99261556, 9.10365258, 8.75129197, 9.34825329, 10.44429514, 11.12300637, 10.65483865, 9.74536512, 9.15844642, 9.47802227, 10.46330973, 11.17151834, 11.12758286, 10.22367137, 9.27834142, 9.30365438, 10.18663464, 11.18485597, 11.15170387, 10.29866589, 9.11982997, 8.91104241, 9.80878568, 10.85327493, 11.20982531, 10.51675446, 9.4423296 , 8.94555971, 9.48191945, 10.33734695, 11.01616257, 10.81084773, 9.62627669, 8.92889768, 9.18181565, 10.05359308, 10.83954893, 10.93661697, 9.91016779, 9.07952059, 9.02910723, 9.85877605, 10.94195393, 11.25544372, 10.71564569, 9.73269725, 9.29290009, 9.89389186, 10.7409496 , 11.16487499, 10.70772978, 9.60656358, 8.96676045, 9.17003758, 10.237347 , 11.07078796, 10.98060372, 10.07576214, 9.16544225, 9.13707365, 9.96672701, 10.93080113, 10.93401823, 10.15791743]

• Your problem I think, but I'm not sure, is that your trend is not really a trend but more of a cycle or a seasonality. I you look at the results of your decomposition, the seasonality component looks more like $f(t) = A(1+sin(t))$ than the trend does. – Skander H. - Reinstate Monica Oct 13 '18 at 19:38
• You are right, my toy example has a seasonality instead of a trend. In my real data, I actually have a mix of both. However, the question is how to fit $\alpha$, regardless of the function $f(t)$. I tried to do it by "removing" $f(t)$ (whether it is a seasonality or trend) but that doesn't seem to work. – RM- Oct 13 '18 at 19:43
• if $\epsilon_t$ are independent of each other, then you can try Ordinary Least Squares; but you would get two coefficient estimates: $\beta$ for $y_t$ and $\alpha$ for $f(t)$. After fitting the model by OLS you could test whether $\beta=1-\alpha$. – javlacalle Oct 16 '18 at 20:34
• If you still want a single coefficient, $\alpha$, you would need to implement the optimisation of the likelihood function by yourself. It shouldn't be that hard, but it may be better to start checking your model using standard procedures already implemented and easier to use. – javlacalle Oct 16 '18 at 20:36
• Does your model comes from a theory or hypothesis that you want to test with real data (maybe $\alpha$ has some meaning in the context of your model)? If not, then you may need to rethink your model or approach. There are specific methods to decompose a time series into trend and seasonal components; state space methods allow defining models with a long-term trend plus a short-term cycle,... Explaining the purpose of your analysis you may get more feedback from others. – javlacalle Oct 17 '18 at 21:00

My suggestion is to parameterise the function $$f(t)$$ and use the methodology of state space models and the Kalman filter. For example, a second order autoregressive, AR(2), process is a relatively general, yet simple, specification that can capture smooth cycles.

Then, you would deal with a Gaussian linear model with an unobserved component. Once you write the state space representation of the model, the Kalman filter can be used to evaluate the likelihood function given the data. By maximising the likelihood function the parameters of the model ($$\alpha$$ and variances $$\sigma^2_\varepsilon$$, $$\sigma^2_\eta$$) can be estimated. In addition, an estimate of the unobserved component $$f(t)$$ is obtained.

Some references about the methodology are given at the end of this answer. Here, I define the model in terms that would allow you to run the Kalman.

The equations of the model:

\begin{align} y_t =& \mu_t + \alpha f_t\\ \mu_t =& (1-\alpha)\mu_{t-1} + \varepsilon_t\,,\quad& \varepsilon_t\sim NID(0,\sigma^2_\varepsilon)\\ f_t =& \phi_1f_{t-1} + \phi_2f_{t-2} + \eta_t\,,\quad& \eta_t\sim NID(0,\sigma^2_\eta)\\ \hbox{Cov}&(\varepsilon,\eta_t)=0 \end{align}

$$y_t$$ is the observed series, which is modeled as of an AR(1) process plus an unobserved transitory component $$f(t)$$. I included two disturbance terms, $$\varepsilon$$ and $$\eta$$, that are independent of each other; the first one can be discarded by setting $$\sigma^2_\varepsilon=0$$.

State space representation (input for the Kalman filter):

$$\begin{eqnarray} y_t=\left( \begin{array}{ccc} 1&\alpha&0 \end{array} \right) \left( \begin{array}{c} \mu_t\\f_t\\f_{t-1} \end{array} \right)\,,\quad \left( \begin{array}{c} \mu_t\\f_t\\f_{t-1} \end{array} \right)= \left( \begin{array}{ccc} 1-\alpha&0&0\\ 0&\phi_1&\phi_2\\ 0&1&0 \end{array} \right) \left( \begin{array}{c} \mu_{t-1}\\f_{t-1}\\f_{t-2} \end{array} \right)\,. \end{eqnarray}$$

The covariance matrix of the disturbance terms is a diagonal matrix (since the disturbances are independent of each other) with diagonal $$(\sigma^2_\varepsilon\,, \sigma^2_\eta\,, 0)$$.

This model is essentially the model proposed by Clark to extract the business cycle from time series of the gross domestic product. I simply added the coefficient $$\alpha$$ in order to resemble your model.

The difference with your original model is that here $$f(t)$$ must be specified by means of a linear model. You may need to try other specifications of $$f(t)$$ that may fit better the overall dynamics of the trend/cycle component in your data. You may also need to include another component to capture, for example, seasonal cycles.

Searching for information about basic structural model you could find a common specification of the seasonal component.

If you cannot come up with a sensible specification that fits your data, then you may need to explore non-parametric methods as I mentioned in a comment.

References:

[1] Clark, P. K. (1987). "The Cyclical Component of U.S. Economic Activity", The Quarterly Journal of Economics, 102, 797-814.

Harvey, A. C. (1989). Forecasting, Structural Time Series Models and the Kalman Filter. Cambridge University Press.

Durbin, J. and Koopman, S. J. (2001). Time Series Analysis by State Space Methods. Oxford University Press.

• This is amazing thanks! About specifying $f_t$, could it be done by measuring trends on the time series? I have several other questions but will read the references first. – RM- Oct 21 '18 at 17:55
• Yes, $f_t$ could be defined as a trend; in fact, an AR(2) process is capable of capturing a trending pattern when $\phi_1+\phi_2$ is close to $1$. You may look for information about the local level model and the Nile data as an example of smooth trends. – javlacalle Oct 23 '18 at 6:56
• Be aware that $\mu$ could also capture a trend, especially when $\alpha$ is close to zero (stochastic trend for $\alpha=0$); some care may be needed in the specification of the model in order to avoid redundant components. – javlacalle Oct 23 '18 at 6:58

From the graphs you have shown in your question, it appears that there may be some periodic signal in your data with a fixed frequency value. To determine whether this is the case, you should generate the Discrete Fourier Transform (DFT) of your data and plot its periodogram. From the plotted periodogram below, it is evident that there is a strong periodic signal in your data, and the estimated frequency is 0.16.

Based on this result, I would recommend using some kind of periodic regression model as your first attempt at modelling the data. Fitting a model with a single sinusoidal wave at this estimated frequency should already explain a lot of the variation in your data. You can add a trend term to this and other terms if you want, but I would start with a simple periodic regression and build it up from there. If you do this, you might find that you get reasonable residuals and you no longer need to worry about any auto-regressive behaviour.

R code: Here is the R code I used to generate these plots:

#Load required libraries
library(ggplot2);
library(stats);

#Import data (put in a text file)
N      <- length(IMPORT$$V1); DATA <- data.frame(Time = 1:N, Value = IMPORT$$V1);

#Create discrete Fourier transform of the data
VALS   <- fft(DATA$$Value - mean(DATA$$Value))[1:(N/2)];
DFT    <- data.frame(Frequency = (0:(N/2-1))/N, Value = VALS, Norm = Mod(VALS));

#Set theme setting for plots
THEME <- list(theme(plot.title    = element_text(hjust = 0.5, face = 'bold', size = 16),
plot.subtitle = element_text(hjust = 0.5, face = 'bold')));

#Generate time series plot
FIGURE1 <- ggplot(data = DATA, aes(x = Time, y = Value)) +
geom_line(size = 1, colour = 'blue') +
THEME +
ggtitle('Plot of time-series data') +
xlab('Time') + ylab('Value');

#Generate periodogram
II <- which(DFT$$Norm == max(DFT$$Norm));
FF <- DFT$Frequency[II]; FIGURE2 <- ggplot(data = DFT, aes(x = Frequency, y = Norm)) + geom_line(size = 1, colour = 'red') + geom_vline(xintercept = FF, size = 1, linetype = 'dashed') + THEME + ggtitle('Periodogram of time-series data') + xlab('Frequency') + ylab('Norm of DFT'); #Plot the time-series and its periodogram FIGURE1; FIGURE2; #Show the estimated frequency of the signal FF; [1] 0.16  • Right, the data was generated by a sin function so there is periodicity. However, the data was only a toy example. In reality I want to fit the$\alpha$parameter instead of looking for alternate models. – RM- Oct 24 '18 at 15:02 • Okay, I didn't realise this wasn't the data set you are using. In any case, your model should conform to your data - it should not be chosen in defiance of the data. The fact that the method I am showing you recovers the actual underlying generation mechanism for the data you presented suggests that it would be an appropriate way to model that data. If your actual data looks different, you should post that instead. – Reinstate Monica Oct 24 '18 at 22:17 • Right. I usually work with unemployment time series data for which there is already a lot of literature and I am sure there are better models. However, I am not so interested in the model but about the questions itself. If there is a time series generated by$y_{t+1} = (1- \alpha) y_t + \alpha f(t) + \epsilon_t$, can we recover$\alpha\$? – RM- Oct 26 '18 at 13:28

This answer is not theoretical but practical and might work in some cases. Use with cautions since it is not guaranteed to work for all cases.

Since $$f(t)$$ change slowly, it is possible to split it into several order 2 polynomials and still have many points. Instead of decomposing the time series in trend + seasonality, we split it and detrend it. After that, we fit the AR(1) model.

We take a rolling window approach and only consider fits that are significant (p-value < 0.001). Then we average over all coefficients. Implementing this methodology I was able to get the following results (in a slightly different example):

real alpha = 0.1

estimated alpha = 0.12591253250291573

standard deviation = 0.08668464697167208

The code I used is provided below.

import pandas as pd
from matplotlib import pylab as plt
from statsmodels.tsa.arima_model import ARMA
import seaborn as sns
import statsmodels.api as sm
import random
import numpy as np
import statsmodels

%matplotlib inline
random.seed(1)

#defining the trend function
def trend(t, amp=1):
return amp*(1 + np.sin(t/10))

#length of time series
n_time_steps = 250
#amplitud of time series
amplitud=10
noise_frac_aplitud= 0.5
#initializing the time series
time_series = np.zeros(n_time_steps)
time_series[0] = trend(0, amplitud)

#The AR(1) parameter. Our goal will be to find this parameter.
alpha = 0.1
#making the time series
for t in range(1,n_time_steps):
time_series[t] = (1 - alpha)*time_series[t - 1] + alpha*trend(t, amp=amplitud) + alpha*np.random.normal(0,noise_frac_aplitud*amplitud)

#passing the time series to a pandas format
dates = sm.tsa.datetools.dates_from_range('2000m1', length=len(time_series))
time_series_pd= pd.Series(time_series, index=dates)

window = 40
n_iter = n_time_steps - window
alpha_list = []
alpha_elite_list = []
#n_time_steps
for i in range(n_iter):
#not rolling window but intervals... to fix
temp_time_series_pd = time_series_pd[i:window + i]
plt.plot(temp_time_series_pd)
res = sm.tsa.detrend(temp_time_series_pd, order=2)
mod = ARMA(res, order=(1,0))
ar1_fit = mod.fit()

score = statsmodels.tsa.arima_model.ARMAResults(mod,ar1_fit.params)

#The alpha parameter is...
alpha_list.append(1 - ar1_fit.params[1])
if score.pvalues[1] < 0.001:
alpha_elite_list.append(1 - ar1_fit.params[1])

print("real alpha = ", alpha)
print("estimated alpha = ", np.mean(alpha_elite_list))
print("standard deviation = ", np.std(alpha_elite_list))