Homework 6
Catherine Dube
2026-03-22
Exercises: 9.1, 9.2, 9.3, 9.5, 9.6, 9.7, 9.8 in Hyndman
Exercise 9.1
Figure 9.32 shows the ACFs for 36 random numbers, 360 random numbers and 1,000 random numbers.
a.
Explain the differences among these figures. Do they all indicate that the data are white noise?
Yes, all of these indicate white noise. The differences between the figures is that as the number of observations goes up, the dashed blue lines representing the confidence interval narrows, and the autocorrelations approach closer to 0.
b
Why are the critical values at different distances from the mean of zero? Why are the autocorrelations different in each figure when they each refer to white noise?
The critical values are at different distances from the mean of zero due to the number of observations. This goes back to Chapter 2 in the book, where these critical values are caluclated by ±1.96/sqrt(T). And T is the number of observations in the time series. As T gets bigger(36 to 360 to 1000), these critical values become narrower.
Exercise 9.2
A classic example of a non-stationary series are stock prices. Plot the daily closing prices for Amazon stock (contained in gafa_stock), along with the ACF and PACF. Explain how each plot shows that the series is non-stationary and should be differenced.
amazon <- gafa_stock %>%
filter(Symbol == 'AMZN')
autoplot(amazon, Close) +
labs(x = "Date", y = "Closing Price", title = "Daily Closing Prices of Amazon Stock")
ACF(amazon, Close) %>%
autoplot() +
labs(title = "ACF for Amazon Stock")## Warning: Provided data has an irregular interval, results should be treated
## with caution. Computing ACF by observation.
PACF(amazon, Close) %>%
autoplot() +
labs(title = "PACF for Amazon Stock")## Warning: Provided data has an irregular interval, results should be treated
## with caution. Computing ACF by observation.
Each plot helps convey that the data is non-stationary and should be differenced.
- The time series plot shows a clear trend in Amazon stock close prices.
- The ACF plot shows a large r1. And from the book I know that for non-stationary data, the value of r1 is often large and positive. Additionally, there are strong autocorrelations that do not drop off quickly (which is they would if this were stationary).
- The PACF plot has a single spike at the beginning.
Exercise 9.3
For the following series, find an appropriate Box-Cox transformation and order of differencing in order to obtain stationary data.
Series
Turkish GDP from global_economy
turkey_gdp <- global_economy %>%
filter(Country == "Turkey")
turkey_gdp %>%
autoplot(GDP) +
labs(title = "Turkish GDP", x = "Year") 
# optimal Box-Cox lambda
lambda_turkey <- turkey_gdp %>%
features(GDP, features = guerrero) %>%
pull(lambda_guerrero)
print(paste0("Lambda: ", lambda_turkey))## [1] "Lambda: 0.157218680108172"# KPSS for # ordinary 1st differences
turkey_gdp %>%
mutate(gdp_trans = box_cox(GDP, lambda_turkey)) %>%
features(gdp_trans, unitroot_ndiffs)## # A tibble: 1 × 2
## Country ndiffs
## <fct> <int>
## 1 Turkey 1# Plot again, transformed and differenced data
turkey_gdp %>%
autoplot(difference(box_cox(GDP, lambda_turkey))) +
labs(title = "Turkish GDP, Transformed and Differenced", x = "Year") ## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
A Box Cox transformation with lambda = 0.1572 and a first differencing will produce stationary data.
Accommodation takings in the state of Tasmania from aus_accommodation
tasmania_accom <- aus_accommodation %>%
filter(State == "Tasmania")
tasmania_accom %>%
autoplot(Takings) +
labs(title = "Accommodation takings in the state of Tasmania", x = "Year") 
# optimal Box-Cox lambda
lambda_tasmania <- tasmania_accom %>%
features(Takings, features = guerrero) %>%
pull(lambda_guerrero)
print(paste0("Lambda: ", lambda_tasmania))## [1] "Lambda: 0.00181964338451332"# Check seasonal differences
tasmania_accom %>%
mutate(takings_trans = box_cox(Takings, lambda_tasmania)) %>%
features(takings_trans, unitroot_nsdiffs) #1## # A tibble: 1 × 2
## State nsdiffs
## <chr> <int>
## 1 Tasmania 1# Check any further ordinary differencing
tasmania_accom %>%
mutate(takings_trans = box_cox(Takings, lambda_tasmania),
takings_sdiff = difference(takings_trans, 4)) %>% # 4 for quarterly data
features(takings_sdiff, unitroot_ndiffs) #none## # A tibble: 1 × 2
## State ndiffs
## <chr> <int>
## 1 Tasmania 0# Plot the transformed and seasonally differenced data
tasmania_accom %>%
autoplot(difference(box_cox(Takings, lambda_tasmania), 4)) +
labs(title = "Transformed and Seasonally Differenced Takings",
y = "Box-Cox(Takings)")## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_line()`).
tasmania_accom %>%
mutate(Trans = box_cox(Takings, lambda_tasmania)) %>%
ACF(difference(Trans, lag=4)) %>%
autoplot() +
labs(title = "ACF Plot: First Seasonal Difference Tasmanian Accomodation Takings")
# Plot the transformed and doubly differenced data (if 1 of each was suggested)
tasmania_accom %>%
autoplot(difference(difference(box_cox(Takings, lambda_tasmania), 4)))+
labs(title = "Transformed, Seasonally Differenced, Ordinarily Differenced Takings",
subtitle = "This is not recommended by the KPSS test, I just wanted to look at it",
y = "Box-Cox(Takings)")## Warning: Removed 5 rows containing missing values or values outside the scale range
## (`geom_line()`).
A Box Cox transformation with lambda = 0.0018 and a first seasonal differencing will produce stationary data. I used the standard methods to get this, but did take a peak at the ACF plot, and it didn’t look definitively stationary because the ACF plot started pretty high but didn’t drop off super quickly.
Monthly sales from souvenirs
# Find the optimal Box-Cox lambda
lambda_souvenirs <- souvenirs %>%
features(Sales, features = guerrero) %>%
pull(lambda_guerrero)
print(paste0("Lambda: ", lambda_souvenirs))## [1] "Lambda: 0.00211822065057588"# Check how many seasonal differences
souvenirs %>%
mutate(sales_trans = box_cox(Sales, lambda_souvenirs)) %>%
features(sales_trans, unitroot_nsdiffs)## # A tibble: 1 × 1
## nsdiffs
## <int>
## 1 1# Check any further ordinary differencing
souvenirs %>%
mutate(sales_trans = box_cox(Sales, lambda_souvenirs),
sales_sdiff = difference(sales_trans, 12)) %>% # 12 for monthly data
features(sales_sdiff, unitroot_ndiffs)## # A tibble: 1 × 1
## ndiffs
## <int>
## 1 0# Plot the transformed and seasonally differenced data
souvenirs %>%
autoplot(difference(box_cox(Sales, lambda_souvenirs), 12)) +
labs(title = "Transformed and Seasonally Differenced Sales",
y = "Box-Cox(Sales)")## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).
souvenirs %>%
mutate(Trans = box_cox(Sales, lambda_souvenirs)) %>%
ACF(difference(Trans, lag=12)) %>%
autoplot() +
labs(title = "ACF Plot: First Seasonal Difference Souvenir Sales")
# Plot the transformed and doubly differenced
souvenirs %>%
autoplot(difference(difference(box_cox(Sales, lambda_souvenirs), 12))) +
labs(title = "Transformed, Seasonally Differenced, Ordinarily Differenced Monthly Sales from Souvenirs",
subtitle = "This is not recommended by the KPSS test, I just wanted to look at it",
y = "Box-Cox(Sales)")## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_line()`).
A Box Cox transformation with lambda = 0.0021 and a first seasonal differencing will produce stationary data. I used the standard methods to get this, but did take a peak at the ACF plot for this example. It didn’t look intuitively stationary to me, but I’m trusting the output from the KPSS test is correct.
Exercise 9.5
For your retail data (from Exercise 7 in Section 2.10), find the appropriate order of differencing (after transformation if necessary) to obtain stationary data.
myseries <- aus_retail %>%
filter(`Series ID` == "A3349476W")
myseries %>% autoplot() + labs(title = "My Series from Exercise 7, Section 2.10")## Plot variable not specified, automatically selected `.vars = Turnover`
# Find the optimal Box-Cox lambda
lambda_my_seriess <- myseries %>%
features(Turnover, features = guerrero) %>%
pull(lambda_guerrero)
print(paste0("Lambda: ", lambda_my_seriess))## [1] "Lambda: 0.113040417623288"# Check how many seasonal differences
myseries %>%
mutate(turnover_trans = box_cox(Turnover, lambda_my_seriess)) %>%
features(turnover_trans, unitroot_nsdiffs)## # A tibble: 1 × 3
## State Industry nsdiffs
## <chr> <chr> <int>
## 1 Queensland Pharmaceutical, cosmetic and toiletry goods retailing 1# Check any further ordinary differencing
myseries %>%
mutate(turnover_trans = box_cox(Turnover, lambda_my_seriess),
turnover_sdiff = difference(turnover_trans, 12)) %>% # 12 for monthly data
features(turnover_sdiff, unitroot_ndiffs)## # A tibble: 1 × 3
## State Industry ndiffs
## <chr> <chr> <int>
## 1 Queensland Pharmaceutical, cosmetic and toiletry goods retailing 0# Plot the transformed and seasonally differenced data
myseries %>%
autoplot(difference(box_cox(Turnover, lambda_my_seriess), 12)) +
labs(title = "Transformed and Seasonally Differenced Sales",
y = "Box-Cox(Sales)")## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).
myseries %>%
mutate(Trans = box_cox(Turnover, lambda_my_seriess)) %>%
ACF(difference(Trans, lag=12)) %>%
autoplot() +
labs(title = "ACF Plot: First Order Difference Souvenirs")
Exercise 9.6
Simulate and plot some data from simple ARIMA models.
a.
Use the following R code to generate data from an AR(1) model with \(\phi_1 = 0.6\) and \(\sigma^2 = 1\). The process starts with \(y_1 = 0\).
y <- numeric(100)
e <- rnorm(100)
for(i in 2:100)
y[i] <- 0.6*y[i-1] + e[i]
sim <- tsibble(idx = seq_len(100), y = y, index = idx)b.
Produce a time plot for the series. How does the plot change as you change $_1 ?
# A few phi values
phi_values <- c(-0.8, -0.5, -0.3, 0, 0.3, 0.5, 0.8)
# Empty list for results
sim_list <- list()
# Iterate though phi values
for (p in phi_values) {
y <- numeric(100)
e <- rnorm(100)
# The AR(1) generation loop
for(i in 2:100) {
y[i] <- p * y[i-1] + e[i]
}
# Store the result as a tibble inside list
sim_list[[as.character(p)]] <- tibble(idx = 1:100, y = y, phi = as.character(p))
}
# Bind into single tsibble
sim_multiple <- bind_rows(sim_list) %>%
as_tsibble(index = idx, key = phi)
# Plots
sim_multiple %>%
autoplot(y) +
facet_wrap(~phi, scales = "free_y") +
labs(title = "AR(1) Simulations for Different Values of Phi",
y = "Value",
x = "Time")
The closer phi is to -1, the more jagged and the more oscillations are crammed together. The closer to 1, the reverse happens.
c.
Write your own code to generate data from an MA(1) model with $_1 = 0.6 and $^2 =1 ?
# Empty y, noise e
y <- numeric(100)
e <- rnorm(100)
# Deal with the first value which just equals its own noise.
y[1] <- e[1]
# Loop
for(i in 2:100) {
y[i] <- e[i] + 0.6 * e[i-1]
}
# tsibble and plot
sim_ma <- tsibble(idx = seq_len(100), y = y, index = idx)d.
Produce a time plot for the series. How does the plot change as you change $_1 ?
# same code as before, swapping phit and theta,AR and MA
# A few the values
theta_values <- c(-0.8, -0.5, -0.3, 0, 0.3, 0.5, 0.8)
# Empty list for results
sim_ma_list <- list()
# Iterate though phi values
for (t in theta_values) {
y <- numeric(100)
e <- rnorm(100)
y[1] <- e[1]
# The MA(1) generation loop
for(i in 2:100) {
y[i] <- e[i] + t * e[i-1]
}
# Store the result as a tibble inside list
sim_ma_list[[as.character(t)]] <- tibble(idx = 1:100, y = y, theta = as.character(t))
}
# Bind into single tsibble
sim_ma_multiple <- bind_rows(sim_ma_list) %>%
as_tsibble(index = idx, key = theta)
# Plots
sim_ma_multiple %>%
autoplot(y) +
facet_wrap(~theta, scales = "free_y") +
labs(title = "MA(1) Simulations for Different Values of Theta",
y = "Value",
x = "Time")
e.
Generate data from an ARMA(1,1) model with $_1 = 0.6, $_1 = 0.6, and $^2 = 1
# Empty y, noise e
y <- numeric(100)
e <- rnorm(100)
# Deal with the first value
y[1] <- e[1]
# ARMA(1,1) Loop
for(i in 2:100) {
# AR part (0.6 * y[i-1]) and MA part (0.6 * e[i-1]) and noise
y[i] <- (0.6 * y[i-1]) + (0.6 * e[i-1]) + e[i]
}
# tsibble
sim_arma <- tsibble(idx = seq_len(100), y = y, index = idx)f.
Generate data from an AR(2) model with $_1 = -0.8, $_2 = 0.3 and $^2 = 1 . . (Note that these parameters will give a non-stationary series.)
# Empty y, noise e
y <- numeric(100)
e <- rnorm(100)
# Deal with the first TWO values bc AR(2) looks back 2 steps
y[1] <- e[1]
y[2] <- -0.8 * y[1] + e[2]
# AR(2) Loop starts at step 3 (because I'm looking back 2 days)
for(i in 3:100) {
y[i] <- -0.8 * y[i-1] + 0.3 * y[i-2] + e[i]
}
# tsibble
sim_ar2 <- tsibble(idx = seq_len(100), y = y, index = idx)g.
Graph the latter two series and compare them.
# tibbles
sim_arma_labeled <- as_tibble(sim_arma) %>% mutate(Model = "ARMA(1,1)")
sim_ar2_labeled <- as_tibble(sim_ar2) %>% mutate(Model = "AR(2)")
# Bind
sim_comparison <- bind_rows(sim_arma_labeled, sim_ar2_labeled) %>%
as_tsibble(index = idx, key = Model)
# Plot
sim_comparison %>%
autoplot(y) +
facet_wrap(~Model, scales = "free_y", ncol = 1) +
labs(title = "Comparison: Stationary ARMA(1,1) vs. Non-Stationary AR(2)",
y = "Value",
x = "Time")
These two are very different. So ARMA(1,1) is a stationary series. It’s hovering about 0 and not straying to any large values. AR(2) is clearly non-stationary. The variance is growing massively over time.
Exercise 9.7
Consider aus_airpassengers, the total number of passengers (in millions) from Australian air carriers for the period 1970-2011.
a.
Use ARIMA() to find an appropriate ARIMA model. What model was selected. Check that the residuals look like white noise. Plot forecasts for the next 10 periods.
# Use ARIMA() to find an appropriate ARIMA model.
fit <- aus_airpassengers %>%
filter(Year < 2011) %>%
model(ARIMA(Passengers))
# What model was selected.
report(fit)## Series: Passengers
## Model: ARIMA(0,2,2)
##
## Coefficients:
## ma1 ma2
## -1.2004 0.3279
## s.e. 0.1922 0.1767
##
## sigma^2 estimated as 4.538: log likelihood=-84.77
## AIC=175.54 AICc=176.23 BIC=180.53# Check that the residuals look like white noise.
fit %>% gg_tsresiduals() +
labs(title = "Residual Diagnostics for Australian Air Passengers")## Warning: `gg_tsresiduals()` was deprecated in feasts 0.4.2.
## ℹ Please use `ggtime::gg_tsresiduals()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# Plot forecasts for the next 10 periods.
fit %>%
forecast(h = 10) %>%
autoplot(aus_airpassengers) +
labs(title = "10-Year Forecast of Australian Air Passengers",
y = "Passengers (Millions)",
x = "Year")
- The selected model was ARIMA(0,2,1).
- The residuals do mostly look like white noise. (The amplitude of them seems to grow slightly over time.) But the ACF shows the autocorrelation well within the confidence bands, and the .resid plot shows the mean at 0, though is appears maybe left skewed slightly.
b.
Write the model in terms of the backshift operator.
The selected model was ARIMA(0,2,1).
- That’s the AR: p=0.
- The I: d=2.
- That’s (1-B)^2 * y_t
- The MA: q=1.
- I know -0.8963 = theta.
- epsilon is the error t
- That’s (1 + thetaB)error or (1-0.8963B)error
The fancy notation:
\[(1-B)^2 y_t = (1 - 0.8963 B) \varepsilon_t\]
c.
Plot forecasts from an ARIMA(0,1,0) model with drift and compare these to part a.
# ARIMA(0,1,0) model with drift
fit_c <- aus_airpassengers %>%
filter(Year < 2011) %>%
model(ARIMA(Passengers ~ 1 + pdq(0,1,0)))
# Check that the residuals look like white noise.
fit_c %>% gg_tsresiduals() +
labs(title = "Residual Diagnostics for ARIMA(0,1,0) with drift")
# Plot new
fit_c %>%
forecast(h = 10) %>%
autoplot(aus_airpassengers) +
labs(title = "Forecast from ARIMA(0,1,0) with drift",
y = "Passengers (Millions)",
x = "Year")
The residuals appear to be less skewed. The ACF shows more positively correlated values than before, whereas it used to appear more equal on both sides. The forecast itself looks to be forecasting slightly lower values too, and is lower than the actual values post 2011.
d.
Plot forecasts from an ARIMA(2,1,2) model with drift and compare these to parts a and c. Remove the constant and see what happens.
# Fit both ARIMA(2,1,2) models (with and without drift)
fit_d <- aus_airpassengers %>%
filter(Year <= 2011) %>%
model(
`ARIMA(2,1,2) w/ drift` = ARIMA(Passengers ~ 1 + pdq(2,1,2)),
`ARIMA(2,1,2) w/o drift` = ARIMA(Passengers ~ 0 + pdq(2,1,2)) # removes the constant
)## Warning: 1 error encountered for ARIMA(2,1,2) w/o drift
## [1] non-stationary AR part from CSSfit_d## # A mable: 1 x 2
## `ARIMA(2,1,2) w/ drift` `ARIMA(2,1,2) w/o drift`
## <model> <model>
## 1 <ARIMA(2,1,2) w/ drift> <NULL model># Plot
fit_d %>%
forecast(h = 10) %>%
autoplot(aus_airpassengers, level = NULL) + # level = NULL hides the confidence cones so the lines are easier to see
labs(title = "Forecasts from ARIMA(2,1,2) With and Without Drift",
y = "Passengers (Millions)",
x = "Year")## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_line()`).
ARIMA(2,1,2) with drift is predicting higher values than the previous two. The one without drift just shows NULL model. This may be because this data has d=1, and it clearly has drift.By forcing a fit without it, the model is supposed to see the difference or growth as 0, which is contradictory to the data, and may make the model coefficients become unstable, resulting in an invalid model, hence the NULL.
e.
Plot forecasts from an ARIMA(0,2,1) model with a constant. What happens?
# Fit the ARIMA(0,2,1) model with a forced constant
fit_e <- aus_airpassengers %>%
filter(Year <= 2011) %>%
model(ARIMA(Passengers ~ 1 + pdq(0,2,1)))## Warning: Model specification induces a quadratic or higher order polynomial trend. This
## is generally discouraged, consider removing the constant or reducing the number
## of differences.# Plot the forecast
fit_e %>%
forecast(h = 10) %>%
autoplot(aus_airpassengers) +
labs(title = "Forecast from ARIMA(0,2,1) with a Constant",
y = "Passengers (Millions)",
x = "Year")
I’m getting a warning from this. “Model specification induces a quadratic or higher order polynomial trend. This is generally discouraged, consider removing the constant or reducing the number of differences.”
Because there are 2 differences, the model accounts for linear trend. If there were only 1 difference, then a constant creates a linear trend. With both of these things, a quadratic equation is happening here, which is generally less suited for forecasting. (The odds are that I don’t actualyl want to forecast something at an accelerating rate like this.)
Exercise 9.8
Consider aus_arrivals, the quarterly number of international visitors to Australia from several countries for the period 1981 Q1 – 2012 Q3.
a.
Select one country and describe the time plot.
# Filter for visitors from the US
us_arrivals <- aus_arrivals %>%
filter(Origin == "US")
# Plot the time series
us_arrivals %>%
autoplot(Arrivals) +
labs(title = "Quarterly Visitors to Australia from the US (1981-2012)",
y = "Number of Arrivals",
x = "Year/Quarter")
The time plot shows an upward trend that loses some velocity around 2000. There are many oscillations that appear to be quarterly seasonality. The amplitude of these oscillations, or variance, appears to be getting bigger over time (though it’s hard to say for certain because they got big in the 1980s as well).
b.
Use differencing to obtain stationary data.
# get num seasonal differences
us_arrivals %>%
features(Arrivals, unitroot_nsdiffs)## # A tibble: 1 × 2
## Origin nsdiffs
## <chr> <int>
## 1 US 1# get num ordinary differences
us_arrivals %>%
mutate(arr_sdiff = difference(Arrivals, 4)) %>% # lag = 4 for quarterly
features(arr_sdiff, unitroot_ndiffs)## # A tibble: 1 × 2
## Origin ndiffs
## <chr> <int>
## 1 US 0# Plot
us_arrivals %>%
autoplot(difference(Arrivals, 4)) +
labs(title = "US Arrivals to Australia (Seasonally Differenced)",
y = "Differenced Arrivals (lag = 4)",
x = "Year/Quarter")## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_line()`).
c.
What can you learn from the ACF graph of the differenced data?
# ACF
us_arrivals %>%
ACF(difference(Arrivals, 4)) %>%
autoplot() +
labs(title = "ACF of Seasonally Differenced US Arrivals",
y = "ACF")
The ACF of the differenced data shows big spikes at the beginning, and then this decays at lag 4. This might mean a non-seasona MA / Moving Average component? (q=1 or q=2). But it’s good that the ACF drops towawrd zero quickly.
d.
What can you learn from the PACF graph of the differenced data?
# ACF
us_arrivals %>%
PACF(difference(Arrivals, 4)) %>%
autoplot() +
labs(title = "PACF of Seasonally Differenced US Arrivals",
y = "PACF")
On the PACF graph, the drop off is still around lag 4, but I see more noticeable negative spikes at the 4 and 8 and somewhat twelve marks (corresponding with seasonality?). However, I also see that this drop off is more gradual, whereas the one in ACF was more sudden. I conclude th is an MA / Moving Average pattern.
e.
What model do these graphs suggest?
ACF has “big spikes at the beginning” and then cuts off; PACF is more of a decay. So I concluded that this suggested MA model. I also know from earlier in the KPSS test that it’s suggested to take 1 seasonal difference and 0 ordinary differences.
So I think the model I would choose would be perhaps ARIMA(0,0,1)(0,1,1) with drift.
f.
Does ARIMA() give the same model that you chose? If not, which model do you think is better?
# ARIMA check
fit_auto <- us_arrivals %>%
model(arima_auto = ARIMA(Arrivals))
# See what it chose
report(fit_auto)## Series: Arrivals
## Model: ARIMA(3,0,3)(0,1,0)[4] w/ drift
##
## Coefficients:
## ar1 ar2 ar3 ma1 ma2 ma3 constant
## -0.2265 -0.4655 0.0961 0.7863 0.9472 0.7735 4414.483
## s.e. 0.1234 0.0925 0.1255 0.0943 0.0551 0.1086 2270.019
##
## sigma^2 estimated as 55384491: log likelihood=-1269.65
## AIC=2555.3 AICc=2556.56 BIC=2577.8ARIMA gave a different model that is more complex.
I suspect that this model is better in terms of accuracy in forecasting (since that that’s the ARIMA() function’s jobs). But it’s less interprettable.
g.
Write the model in terms of the backshift operator, then without using the backshift operator.
Trying to break down the components:
- AR part (\(p=3\)): \((1 - \phi_1 B - \phi_2 B^2 - \phi_3 B^3)\)
- Seasonal Difference (\(D=1\)): \((1 - B^4)\)
- MA part (\(q=3\)): \((1 + \theta_1 B + \theta_2 B^2 + \theta_3 B^3)\)
- Constant: \(c\) (this is 4414.483 from the model)
- Error: \(\varepsilon_t\)$
The Equation (Backshift Notation):
\[(1 - \phi_1 B - \phi_2 B^2 - \phi_3 B^3)(1 - B^4) y_t = c + (1 + \theta_1 B + \theta_2 B^2 + \theta_3 B^3) \varepsilon_t\]
The Equation (Without Backshift):
To write it without \(B\), I’ll distribute the terms. The workflow is below (jump to the end to get the final equation because I just could not do it in my head):
Backshift Equation:
\[(1 - \phi_1 B - \phi_2 B^2 - \phi_3 B^3)(1 - B^4) y_t = c + (1 + \theta_1 B + \theta_2 B^2 + \theta_3 B^3) \varepsilon_t\]
Deal with the Right Hand Side (MA):
Distribute the error term \(\varepsilon_t\) through the MA polynomial. The backshift operator \(B^k\) shifts the error back \(k\) periods: \[c + \varepsilon_t + \theta_1 B \varepsilon_t + \theta_2 B^2 \varepsilon_t + \theta_3 B^3 \varepsilon_t\] \[= c + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \theta_2 \varepsilon_{t-2} + \theta_3 \varepsilon_{t-3}\]
Deal with the Left hand side (Differenicng + AR):
Applying seasonal difference \((1 - B^4)\) directly to \(y_t\): \[(1 - B^4) y_t = y_t - y_{t-4}\]
Substituting this into the AR polynomial: \[(1 - \phi_1 B - \phi_2 B^2 - \phi_3 B^3)(y_t - y_{t-4})\]
Distributing \((y_t - y_{t-4})\) through the AR terms: \[1(y_t - y_{t-4}) - \phi_1 B(y_t - y_{t-4}) - \phi_2 B^2(y_t - y_{t-4}) - \phi_3 B^3(y_t - y_{t-4})\]
Applying the backshift operator \(B\) to the terms inside the parentheses (shifting them back 1, 2, and 3 periods respectively): \[(y_t - y_{t-4}) - \phi_1(y_{t-1} - y_{t-5}) - \phi_2(y_{t-2} - y_{t-6}) - \phi_3(y_{t-3} - y_{t-7})\]
Solve for \(y_t\):
Setting the expanded Left Hand Side equal to the expanded Right Hand Side: \[(y_t - y_{t-4}) - \phi_1(y_{t-1} - y_{t-5}) - \phi_2(y_{t-2} - y_{t-6}) - \phi_3(y_{t-3} - y_{t-7}) = c + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \theta_2 \varepsilon_{t-2} + \theta_3 \varepsilon_{t-3}\]
Moving other terms over to solve for \(y_t\)
Final Equation: \[y_t = y_{t-4} + \phi_1(y_{t-1} - y_{t-5}) + \phi_2(y_{t-2} - y_{t-6}) + \phi_3(y_{t-3} - y_{t-7}) + c + \varepsilon_t + \theta_1\varepsilon_{t-1} + \theta_2\varepsilon_{t-2} + \theta_3\varepsilon_{t-3}\]