The time series \(\{x_t\}\) is an autoregressive process of order \(p\), denoted as \(AR(p)\), if \[
x_t = \alpha_1 x_{t-1} + \alpha_2 x_{t-2} + \alpha_3 x_{t-3} + \cdots + \alpha_{p-1} x_{t-(p-1)} + \alpha_p x_{t-p} + w_t ~~~~~~~~~~~~~~~~~~~~~~~ (4.15)
\]
where \(\{w_t\}\) is white noise and the \(\alpha_i\) are the model parameters with \(\alpha_p \ne 0\).
The \(AR(p)\) model can be expressed as: \[
\underbrace{\left( 1 - \alpha_1 \mathbf{B} - \alpha_2 \mathbf{B}^2 - \cdots - \alpha_p \mathbf{B}^p \right)}_{\theta_p \left( \mathbf{B} \right)} x_t = w_t
\]
Moving Average (MA) Models
The definition of an \(MA(q)\) model is:
Definition of a Moving Average (MA) Model
We say that a time series \(\{x_t\}\) is a moving average process of order \(q\), denoted as \(MA(q)\), if each term in the time series is a linear combination of the current white noise term and the \(q\) most recent past white noise terms.
It is given as: \[
x_t = w_t + \beta_1 w_{t-1} + \beta_2 w_{t-2} + \beta_3 w_{t-3} + \cdots + \beta_{q-1} w_{t-(q-1)} + \beta_q w_{t-q}
\]
where \(\{w_t\}\) is white noise with zero mean and variance \(\sigma_w^2\), and the \(\beta_i\) are the model parameters with \(\beta_q \ne 0\).
Written in terms of the backward shift operator, we have \[
x_t = \underbrace{\left( 1 + \beta_1 \mathbf{B} + \beta_2 \mathbf{B}^2 + \beta_3 \mathbf{B}^3 + \cdots + \beta_{q-1} \mathbf{B}^{q-1} + \beta_q \mathbf{B}^{q} \right)}_{\phi_q(\mathbf{B})} w_t
\]
Putting the \(AR\) and \(MA\) models together, we get the \(ARMA\) model.
Definition of an Autogregressive Moving Average (ARMA) Model
A time series \(\{ x_t \}\) follows an autoregressive moving average (ARMA) model of order \((p, q)\), which we write as \(ARMA(p,q)\), if it can be written as:
We can write this as: \[
\theta_p \left( \mathbf{B} \right) x_t
=
\phi_q \left( \mathbf{B} \right) w_t
\]
Facts about ARMA Processes
The following facts are true for \(ARMA(p,q)\) processes:
The ARMA process is stationary if all the roots of \(\theta_p \left( \mathbf{B} \right)\) are greater than 1 in absolute value.
The ARMA process is invertible if all the roots of \(\phi_q \left( \mathbf{B} \right)\) are greater than 1 in absolute value.
The special case \(ARMA(p,0)\) is the \(AR(p)\) model.
The special case \(ARMA(0,q)\) is the \(MA(q)\) model.
An \(ARMA\) model will usually require fewer parameters than a single \(MA\) or \(AR\) model. This is called parameter parsimony.
If \(\theta\) and \(\phi\) have a common factor, a stationary model can be simplified. This is called parameter redundancy. As an example, the model \[
\left( 1 - \frac{1}{2} \mathbf{B} \right)\left( 1 - \frac{1}{3} \mathbf{B} \right) x_t
=
\left( 1-\frac{1}{2} \mathbf{B} \right)\left( 1 - \frac{1}{4} \mathbf{B} \right) w_t
\] is the same as the model \[
\left( 1 - \frac{1}{3} \mathbf{B} \right) x_t
=
\left( 1 - \frac{1}{4} \mathbf{B} \right) w_t
\]
Comparison of AR and MA Models
ACF and PACF of an \(AR(p)\) Process
We can use the pacf and acf plots to assess if an \(AR(p)\) or \(MA(q)\) model is appropriate. For an \(AR(p)\) or \(MA(q)\) process, we observe the following:
AR(p)
MA(q)
ARMA(p,q)
ACF
Tails off
Cuts off after lag \(q\)
Tails off
PACF
Cuts off after lag \(p\)
Tails off
Tails off
Analyzing Time Series with a Regular Seasonal Pattern
Here are some steps you can use to guide your work as you analyze time series with a regular seasonal pattern. Even though these are presented linearly, you may find it helpful to iterate between some of the steps as needed.
Analyzing a Time Series with a Regular Seasonal Pattern
Create a time plot of the series
Determine an appropriate model for the trend and seasonality; call this Model 1
Regression
Holt-Winters
Additive / Multiplicative Decomposition
Other Techniques
Use AIC/AICc/BIC (or other metrics) and your reasoning skills to choose the best model
Obtain the residuals
Check the residuals of Model 1 for evidence of autocorrelation and non-stationarity
Time plot of the residuals
ACF plot
PACF plot
Fit an AR/MA/ARMA model to the (stationary) residuals; call this Model 2
Check the residuals of Model 2 for evidence of autocorrelation and normality
ACF plot
PACF plot
Histogram
Forecast the trend and seasonality of the time series using Model 1
Forecast the residuals of Model 1 using Model 2
Add the two forecasts together
Plot and summarize the resulting forecast
Refine the models above as needed
Class Activity: Model for the Residuals from the Rexburg Weather Model (15 min)
Show the code
weather_df <- rio::import("https://byuistats.github.io/timeseries/data/rexburg_weather_monthly.csv") |>mutate(dates =my(date_text)) |>filter(dates >=my("1/2008") & dates <=my("12/2023")) |>rename(x = avg_daily_high_temp) |>mutate(TIME =1:n()) |>mutate(zTIME = (TIME -mean(TIME)) /sd(TIME),month_index =month(dates),month_index =factor(month_index, levels =1:12, labels = month.abb)) |>as_tsibble(index = TIME)weather_df |>as_tsibble(index = dates) |>autoplot(.vars = x) +geom_smooth(method ="lm", se =FALSE, color ="#F0E442") +labs(x ="Month",y ="Mean Daily High Temperature (Fahrenheit)",title ="Time Plot of Mean Daily Rexburg High Temperature by Month",subtitle =paste0("(", format(weather_df$dates %>%head(1), "%b %Y"), endash, format(weather_df$dates %>%tail(1), "%b %Y"), ")") ) +theme_minimal() +theme(plot.title =element_text(hjust =0.5),plot.subtitle =element_text(hjust =0.5) )
Figure 1: Time plot of the monthly mean daily high temperatures in Rexburg Idaho, in degrees Fahrenheit
However \(\beta_1\) is for our normalized time variable. If we want to put it in terms of changes each month we would divide our estimate by \(\sigma_t\) (getting an estimate of 0.0124). If we would want it in terms of changes each year we would multiply by 12 (getting an estimate of 0.148).
Figure 2: Time plot of the monthly mean daily high temperatures in Rexburg Idaho, in degrees Fahrenheit; the fitted values from the regression model are given in orange
Fitting an ARMA(p,q) Model
First, we create an acf and pacf plot of the residuals from the model above.
Figure 28: Partial correlogram of the residuals from the fitted regression model
Show the code
################################################################# Be sure to adjust this model to match what you did above ################################################################### Model without any seasonal trendsretail_no_seasonal_lm <- retail_ts |>model(retail_full_quad =TSLM(log(sales_millions) ~ (in_great_recession * std_t) + (after_great_recession * std_t) + (in_great_recession *I(std_t^2) ) + (after_great_recession *I(std_t^2) )+ (in_great_recession *I(std_t^3) ) + (after_great_recession *I(std_t^3) ) ) )retail_ts |>autoplot(.vars = sales_millions) +# ggplot(mapping = aes(x = month, y = sales_millions)) +# geom_line() +geom_line(data =augment(retail_no_seasonal_lm),aes(x = month, y = .fitted),color ="#E69F00",linewidth =1 ) +geom_line(data =augment(retail_final_lm),aes(x = month, y = .fitted),color ="blue",alpha =0.7,# linewidth = 1 ) +labs(x ="Month",y ="Sales (Millions of U.S. Dollars)",title =paste0(retail_ts$business[1], " (", retail_ts$naics[1], ")") ) +theme_minimal() +theme(plot.title =element_text(hjust =0.5))
Figure 29: Time plot of the retail sales data for full-service restaurants; the regression curve without the seasonal terms is given in orange; the predicted values from the model are given in blue
Table 3: Results from the process of fitting an ARMA model to the residuals of the model for the total sales in full-service restaurants in the United States
