Homework 6
Seung Min Song
2024-03-17
9.1, 9.2, 9.3, 9.5, 9.6, 9.7, 9.8 in Hyndman
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? Left: ACF for a white noise series of 36 numbers. Middle: ACF for a white noise series of 360 numbers. Right: ACF for a white noise series of 1,000 numbers. Figure 9.32: Left: ACF for a white noise series of 36 numbers. Middle: ACF for a white noise series of 360 numbers. Right: ACF for a white noise series of 1,000 numbers.
A: On the left, some bars are slightly protruding from the boundary, which may suggest some autocorrelation. However, this could also be due to chance, given the small sample size.
In the middle plot, as the sample size increases, the autocorrelation values at various lags converge towards zero, and more of them fall within the significance range. This pattern better reflects a white noise process compared to the leftmost plot.
In the rightmost graph, with large sample sizes, the autocorrelation values are very close to zero and fall within the significance range for all lags, indicating a strong presence of a white noise process.
All three figures aim to illustrate characteristics of white noise data.
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?
A:
Thresholds vary in their distance from mean zero, primarily due to the sample size of the time series data. Despite all being referenced to white noise, differences in autocorrelation for each figure may be due to random variation and the effect of sample size.
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_stock <- gafa_stock %>%
filter(Symbol == "AMZN")
head(amazon_stock)## # A tsibble: 6 x 8 [!]
## # Key: Symbol [1]
## Symbol Date Open High Low Close Adj_Close Volume
## <chr> <date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 AMZN 2014-01-02 399. 399. 394. 398. 398. 2137800
## 2 AMZN 2014-01-03 398. 403. 396. 396. 396. 2210200
## 3 AMZN 2014-01-06 396. 397 388. 394. 394. 3170600
## 4 AMZN 2014-01-07 395. 398. 394. 398. 398. 1916000
## 5 AMZN 2014-01-08 398. 403 396. 402. 402. 2316500
## 6 AMZN 2014-01-09 404. 407. 398. 401. 401. 2103000# Plotting the daily closing prices for Amazon
amazon_stock %>%
ggplot(aes(x = Date, y = Close)) +
geom_line() +
ggtitle("Daily Closing Prices for Amazon") +
xlab("Date") +
ylab("Closing Price")
par(mfrow = c(2, 1))
# Plotting the ACF for Amazon's closing stock prices
acf(amazon_stock$Close, main="ACF for Amazon Closing Prices")
# Plotting the PACF for Amazon's closing stock prices
pacf(amazon_stock$Close, main="PACF for Amazon Closing Prices")
If the mean changes with time, it is a nonstationary time series. The ACF plot shows a very gradual decrease in the autocorrelation coefficient as the delay increases. For stationary time series, this correlation is expected to quickly fall to zero as lag increases. Instead, this slow decline suggests that the autocorrelation persists over time, which is a common feature of nonstationary series.
PACF Plot for Amazon Closing Prices: The PACF plot does not show a clear cut-off point after which the partial autocorrelations are insignificant (within the blue dashed confidence bounds). If the series were stationary, we would typically see a significant spike at small lags (for example, at lag 1 or 2) and then the PACF would fall within the confidence interval for the remaining lags, indicating no significant partial correlation. However, the PACF here does not provide as clear evidence of non-stationarity as the ACF does.
After differencing, the ACF and PACF plots would be re-examined. For a stationary differenced series, you would expect the ACF to drop to insignificance more quickly, and the PACF might show a significant spike at the first lag if the data follow an AR(1) process after differencing. If the data are overdifferenced, you might see a significant negative spike at the first lag in the ACF.
9.3
For the following series, find an appropriate Box-Cox transformation and order of differencing in order to obtain stationary data.
a.
Turkish GDP from global_economy.
turkey_gdp <- global_economy %>%
filter(Country == "Turkey") %>%
select(Year, GDP)
lambda_opt <- BoxCox.lambda(turkey_gdp$GDP)
turkey_gdp_transformed <- BoxCox(turkey_gdp$GDP, lambda_opt)
turkey_gdp_diff <- diff(turkey_gdp_transformed)
par(mfrow=c(2,2))
plot(turkey_gdp$Year, turkey_gdp$GDP, type="l",
main="Non-transformed Turkish GDP")
plot(turkey_gdp$Year[-1], turkey_gdp_diff, type="l",
main=paste("Differenced Turkish GDP with λ =", round(lambda_opt, 2)))
acf(turkey_gdp$GDP, main="ACF for Non-transformed Turkish GDP")
pacf(turkey_gdp$GDP, main="PACF for Non-transformed Turkish GDP")
acf(turkey_gdp_diff, main=paste("ACF for Differenced Turkish GDP\n with λ =", round(lambda_opt, 2)))
pacf(turkey_gdp_diff, main=paste("PACF for Differenced Turkish GDP\n with λ =", round(lambda_opt, 2)))
b.
Accommodation takings in the state of Tasmania from aus_accommodation.
Tasmania <- aus_accommodation %>%
filter(State == "Tasmania") %>%
select(Date, Takings)
lambda_opt <- BoxCox.lambda(Tasmania$Takings)
lambda_opt## [1] -0.005076712Tasmania_transformed <- BoxCox(Tasmania$Takings, lambda_opt)
Tasmania_diff <- diff(Tasmania_transformed)par(mfrow=c(2,2))
plot(Tasmania$Date, Tasmania$Takings, type="l",
main="Non-transformed Tasmania Takings")
plot(Tasmania$Date[-1], Tasmania_diff, type="l",
main=paste("Differenced Tasmania Takings\n with λ =", round(lambda_opt, 2)))
acf(Tasmania$Takings, main="ACF for Non-transformed Takings")
pacf(Tasmania$Takings, main="PACF for Non-transformed Takings")
acf(Tasmania_diff, main=paste("ACF for Differenced Tasmania Takings\n with λ =", round(lambda_opt, 2)))
pacf(Tasmania_diff, main=paste("PACF for Differenced Tasmania Takings\n with λ =", round(lambda_opt, 2)))
C.
Monthly sales from souvenirs.
souvenirs## # A tsibble: 84 x 2 [1M]
## Month Sales
## <mth> <dbl>
## 1 1987 Jan 1665.
## 2 1987 Feb 2398.
## 3 1987 Mar 2841.
## 4 1987 Apr 3547.
## 5 1987 May 3753.
## 6 1987 Jun 3715.
## 7 1987 Jul 4350.
## 8 1987 Aug 3566.
## 9 1987 Sep 5022.
## 10 1987 Oct 6423.
## # ℹ 74 more rowslambda_opt <- BoxCox.lambda(souvenirs$Sales)
lambda_opt## [1] -0.2444328souvenirs_transformed <- BoxCox(souvenirs$Sales, lambda_opt)
souvenirs_diff <- diff(souvenirs_transformed)
par(mfrow=c(2,2))
plot(souvenirs$Month, souvenirs$Sales, type="l",
main="Non-transformed Souvenirs Sales")
plot(souvenirs$Month[-1], souvenirs_diff, type="l",
main=paste("Differenced Souvenirs Sales\n with λ =", round(lambda_opt, 2)))
acf(souvenirs$Sales, main="ACF for Non-transformed Souvenirs\n Sales")
pacf(souvenirs$Sales, main="PACF for Non-transformed Souvenirs\n Sales")
acf(souvenirs_diff, main=paste("ACF for Differenced Souvenirs Sales\n with λ =", round(lambda_opt, 2)))
pacf(souvenirs_diff, main=paste("PACF for Differenced Souvenirs Sales\n with λ =", round(lambda_opt, 2)))
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.
set.seed(1234)
myseries <- aus_retail %>%
filter(`Series ID` == sample(aus_retail$`Series ID`,1))
myseries## # A tsibble: 441 x 5 [1M]
## # Key: State, Industry [1]
## State Industry `Series ID` Month Turnover
## <chr> <chr> <chr> <mth> <dbl>
## 1 Tasmania Cafes, restaurants and takeaway food … A3349520V 1982 Apr 5.4
## 2 Tasmania Cafes, restaurants and takeaway food … A3349520V 1982 May 5.5
## 3 Tasmania Cafes, restaurants and takeaway food … A3349520V 1982 Jun 5.1
## 4 Tasmania Cafes, restaurants and takeaway food … A3349520V 1982 Jul 5.5
## 5 Tasmania Cafes, restaurants and takeaway food … A3349520V 1982 Aug 5.5
## 6 Tasmania Cafes, restaurants and takeaway food … A3349520V 1982 Sep 5.7
## 7 Tasmania Cafes, restaurants and takeaway food … A3349520V 1982 Oct 5.9
## 8 Tasmania Cafes, restaurants and takeaway food … A3349520V 1982 Nov 5.9
## 9 Tasmania Cafes, restaurants and takeaway food … A3349520V 1982 Dec 6.7
## 10 Tasmania Cafes, restaurants and takeaway food … A3349520V 1983 Jan 5.5
## # ℹ 431 more rowsUse kpass to determine the appropriate differential order for the Turnover column. Looking at the results, the number was 1.
ndiffs_diff <- ndiffs(myseries$Turnover, test = "kpss")
ndiffs_diff## [1] 1myseries_diff <- diff(myseries$Turnover, differences = ndiffs_diff)par(mfrow=c(2,2))
plot(myseries$Month, myseries$Turnover, type="l",
main="Non-transformed Turnover")
plot(myseries_diff, main=paste("Differenced Data with λ =", round(lambda_opt, 2)), xlab="Time", ylab="Transformed and Differenced Turnover")
acf(myseries_diff, main=paste("ACF of Differenced Data with λ =", round(lambda_opt, 2)))
pacf(myseries_diff, main=paste("PACF of Differenced Data with λ =", round(lambda_opt, 2)))
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 ϕ1=0.6 and σ2=1. The process starts with y1=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 ?
These values are chosen because they are within the stationary range for an AR(1) model, which is −1<ϕ1<1. -1 or closer to 1 indicates more extreme behavior, but is still within the range where the AR(1) process stalls.
simulate_ar1 <- function(phi, n = 100, sigma = 1) {
y <- numeric(n)
e <- rnorm(n, mean = 0, sd = sigma)
y[1] <- rnorm(1, mean = 0, sd = sigma / sqrt(1 - phi^2))
for(i in 2:n) {
y[i] <- phi * y[i - 1] + e[i]
}
return(y)
}
set.seed(1234)
par(mfrow = c(2, 2))
phi_values <- c(0.5, -0.5, 0.9, -0.9)
for (phi in phi_values) {
y <- simulate_ar1(phi)
plot(y, type = "l", main = paste("AR(1) Series with phi =", phi),
xlab = "Time", ylab = "Value")
}
c.
Write your own code to generate data from an MA(1) model with θ1=0.6 and σ2=1 .
In fact, the first item in the e array is unused, so we don’t compute a value for e[0], producing a total of 101 error terms from e[1] to e[101]. So we can run the for loop until 2:101 to calculate all 100 values in the y array.
set.seed(123)
theta1 <- 0.6
sigma2 <- 1
e <- rnorm(101, sd = sqrt(sigma2))
y <- numeric(100)
for (i in 2:101) {
y[i-1] <- e[i] + theta1 * e[i-1]
}
plot(y, type="o", main="Simulated MA(1) Data with θ1 = 0.6", ylab="Value", xlab="Time")
d.
Produce a time plot for the series. How does the plot change as you change θ1 ?
The value and sign of θ1 have a significant impact on the behavior of the MA(1) process. The closer the θ1 value is to 1 or -1, the stronger the persistence or reversal.
set.seed(123)
generate_MA1 <- function(theta1, n = 100, sigma2 = 1) {
e <- rnorm(n + 1, mean = 0, sd = sqrt(sigma2))
y <- numeric(n)
for(i in 1:n) {
y[i] <- e[i + 1] + theta1 * e[i]
}
return(y)
}
theta1_values <- c(0.6, -0.6, 0.9, -0.9)
par(mfrow=c(2,2))
for(theta1 in theta1_values) {
y <- generate_MA1(theta1)
plot(y, type="l", main=paste("MA(1) Series with θ1 =", theta1), xlab="Time", ylab="Value")
}
e.
Generate data from an ARMA(1,1) model with ϕ1= , θ1=0.6 and σ2=1 .
set.seed(123)
phi1_arma <- 0.6
theta1_arma <- 0.6
sigma2_arma <- 1
e_arma <- rnorm(102, mean = 0, sd = sqrt(sigma2_arma))
y_arma <- numeric(102)
y_arma[1] <- e_arma[1] / sqrt(1 - phi1_arma^2)
for (i in 3:102) {
y_arma[i] <- phi1_arma * y_arma[i - 1] + e_arma[i] + theta1_arma * e_arma[i - 1]
}
y_arma <- y_arma[-c(1:2)] 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.)
set.seed(123)
phi1_ar2 <- -0.8
phi2_ar2 <- 0.3
sigma2_ar2 <- 1
e_ar2 <- rnorm(102, mean = 0, sd = sqrt(sigma2_ar2))
y_ar2 <- numeric(102)
y_ar2[1] <- e_ar2[1]
y_ar2[2] <- e_ar2[2]
for (i in 3:102) {
y_ar2[i] <- phi1_ar2 * y_ar2[i - 1] + phi2_ar2 * y_ar2[i - 2] + e_ar2[i]
}
y_ar2 <- y_ar2[-c(1:2)] g.
Graph the latter two series and compare them.
- The ARMA(1,1) model represents a stable process, while the AR(2) model tends to be unstable.
- The ARMA(1,1) model shows a process that stably oscillates around the average, while the AR(2) model shows a tendency to oscillate from the beginning and then diverge rapidly.
With the given parameters, it shows that the ARMA(1,1) model is suitable for modeling stable time series, and the AR(2) model is not suitable for modeling unstable time series.
par(mfrow=c(2, 1))
plot(y_arma, type="l", main="Simulated ARMA(1,1) Data with phi1 = 0.6 and theta1 = 0.6",
ylab="Value", xlab="Time")
plot(y_ar2, type="l", main="Simulated AR(2) Data with phi1 = -0.8 and phi2 = 0.3",
ylab="Value", xlab="Time")
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.
The ACF plot shows the autocorrelation function (ACF) of the residuals. There is no autocorrelation outside the dashed blue line that represents the confidence interval. The model’s residuals typically appear to behave like white noise, indicating that the model is a good fit to the data.
aus_airpassengers## # A tsibble: 47 x 2 [1Y]
## Year Passengers
## <dbl> <dbl>
## 1 1970 7.32
## 2 1971 7.33
## 3 1972 7.80
## 4 1973 9.38
## 5 1974 10.7
## 6 1975 11.1
## 7 1976 10.9
## 8 1977 11.3
## 9 1978 12.1
## 10 1979 13.0
## # ℹ 37 more rowsaus_passengers_ts <- ts(aus_airpassengers$Passengers, start = min(aus_airpassengers$Year), frequency = 1)
fit <- auto.arima(aus_passengers_ts)
summary(fit)## Series: aus_passengers_ts
## ARIMA(0,2,1)
##
## Coefficients:
## ma1
## -0.8963
## s.e. 0.0594
##
## sigma^2 = 4.308: log likelihood = -97.02
## AIC=198.04 AICc=198.32 BIC=201.65
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.3453292 2.008183 1.215008 0.9528709 4.809966 0.699953 -0.1743841checkresiduals(fit)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,2,1)
## Q* = 6.5774, df = 8, p-value = 0.5828
##
## Model df: 1. Total lags used: 9forecast_fit <- forecast(fit, h = 10)
plot(forecast_fit)
b.
Write the model in terms of the backshift operator.
(1−B)2 Yt =(1+θ1B)ϵt
(1−B) 2 Y t =(1−0.8963B)ϵ tc.
Plot forecasts from an ARIMA(0,1,0) model with drift and compare these to part a.
drift is used because it is constant trend over time in the data. Both graph’s drift are quite similar. They both show an upward trend, which is expected since the data likely exhibits a long-term increase over time.
fit_drift <- Arima(aus_passengers_ts, order=c(0,1,0), include.drift=TRUE)
summary(fit_drift)## Series: aus_passengers_ts
## ARIMA(0,1,0) with drift
##
## Coefficients:
## drift
## 1.4191
## s.e. 0.3014
##
## sigma^2 = 4.271: log likelihood = -98.16
## AIC=200.31 AICc=200.59 BIC=203.97
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.0001255232 2.022083 1.30933 -2.459026 5.78925 0.7542905
## ACF1
## Training set -0.02389837checkresiduals(fit_drift)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,0) with drift
## Q* = 6.7444, df = 9, p-value = 0.6637
##
## Model df: 0. Total lags used: 9forecast_fit_drift <- forecast(fit_drift, h = 10)
plot(forecast_fit_drift)
par(mfrow = c(2, 1))
plot(forecast_fit, main="Forecast from auto.arima")
plot(forecast_fit_drift, main="Forecast from ARIMA(0,1,0) with drift")
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.
arima_212_with_drift <- Arima(aus_passengers_ts, order=c(2,1,2), include.drift=TRUE)
summary(arima_212_with_drift)## Series: aus_passengers_ts
## ARIMA(2,1,2) with drift
##
## Coefficients:
## ar1 ar2 ma1 ma2 drift
## -0.5518 -0.7327 0.5895 1.0000 1.4101
## s.e. 0.1684 0.1224 0.0916 0.1008 0.3162
##
## sigma^2 = 4.031: log likelihood = -96.23
## AIC=204.46 AICc=206.61 BIC=215.43
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.006116278 1.875262 1.181744 -2.136362 5.276224 0.6807896
## ACF1
## Training set -0.02039146checkresiduals(arima_212_with_drift)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,1,2) with drift
## Q* = 4.9999, df = 5, p-value = 0.4159
##
## Model df: 4. Total lags used: 9forecast_arima_212_with_drift <- forecast(arima_212_with_drift, h = 10)
arima_212_no_constant <- Arima(aus_passengers_ts, order=c(2,1,2), seasonal=c(0,0,0), include.drift=FALSE, method="ML")
summary(arima_212_no_constant )## Series: aus_passengers_ts
## ARIMA(2,1,2)
##
## Coefficients:
## ar1 ar2 ma1 ma2
## 1.8081 -0.8111 -1.9642 1.0000
## s.e. 0.1121 0.1126 0.1091 0.1102
##
## sigma^2 = 3.9: log likelihood = -97.19
## AIC=204.38 AICc=205.88 BIC=213.52
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.4037228 1.866736 1.196676 1.255504 4.692374 0.6893919
## ACF1
## Training set -0.06106455checkresiduals(arima_212_no_constant)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,1,2)
## Q* = 4.2295, df = 5, p-value = 0.5169
##
## Model df: 4. Total lags used: 9forecast_arima_212_no_constant <- forecast(arima_212_no_constant, h = 10)
par(mfrow=c(2,1))
plot(forecast_fit, main="Forecast from auto.arima")
plot(forecast_fit_drift, main="Forecast from ARIMA(0,1,0) with drift")
plot(forecast_arima_212_with_drift, main="Forecast from ARIMA(2,1,2) with drift")
plot(forecast_arima_212_no_constant, main="Forecast from ARIMA(2,1,2) without constant")
e.
Plot forecasts from an ARIMA(0,2,1) model with a constant. What happens?
ARIMA(0,2,1) is essentially a second-difference model with a moving average component.If there is a clear trend in the time series adding a constant to a model that already has differences may not significantly change the model’s fit to the data. The mean change is already small or negligible after differencing.
arima_with_constant <- Arima(aus_passengers_ts, order=c(0,2,1), include.constant=TRUE)
summary(arima_with_constant)## Series: aus_passengers_ts
## ARIMA(0,2,1)
##
## Coefficients:
## ma1
## -0.8963
## s.e. 0.0594
##
## sigma^2 = 4.308: log likelihood = -97.02
## AIC=198.04 AICc=198.32 BIC=201.65
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.3453292 2.008183 1.215008 0.9528709 4.809966 0.699953 -0.1743841arima_without_constant <- Arima(aus_passengers_ts, order=c(0,2,1), include.constant=FALSE)
summary(arima_without_constant)## Series: aus_passengers_ts
## ARIMA(0,2,1)
##
## Coefficients:
## ma1
## -0.8963
## s.e. 0.0594
##
## sigma^2 = 4.308: log likelihood = -97.02
## AIC=198.04 AICc=198.32 BIC=201.65
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.3453292 2.008183 1.215008 0.9528709 4.809966 0.699953 -0.17438419.8
For the United States GDP series (from global_economy):
a.
if necessary, find a suitable Box-Cox transformation for the data.
An increasing spread means that the change is proportional to the level of GDP. The Box-Cox transformation can stabilize the variance.
US_economy <- global_economy %>%
filter(Code == "USA")
US_economy## # A tsibble: 58 x 9 [1Y]
## # Key: Country [1]
## Country Code Year GDP Growth CPI Imports Exports Population
## <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 United States USA 1960 5.43e11 NA 13.6 4.20 4.97 180671000
## 2 United States USA 1961 5.63e11 2.30 13.7 4.03 4.90 183691000
## 3 United States USA 1962 6.05e11 6.10 13.9 4.13 4.81 186538000
## 4 United States USA 1963 6.39e11 4.40 14.0 4.09 4.87 189242000
## 5 United States USA 1964 6.86e11 5.80 14.2 4.10 5.10 191889000
## 6 United States USA 1965 7.44e11 6.40 14.4 4.24 4.99 194303000
## 7 United States USA 1966 8.15e11 6.50 14.9 4.55 5.02 196560000
## 8 United States USA 1967 8.62e11 2.50 15.3 4.63 5.05 198712000
## 9 United States USA 1968 9.42e11 4.80 16.0 4.94 5.08 200706000
## 10 United States USA 1969 1.02e12 3.10 16.8 4.95 5.09 202677000
## # ℹ 48 more rowsplot(US_economy$Year, US_economy$GDP, type="l",
main="US GDP")
lambda_opt <- BoxCox.lambda(US_economy$GDP)
US_economy$GDP_transformed <- BoxCox(US_economy$GDP, lambda_opt)
plot(US_economy$Year, US_economy$GDP_transformed, type="l", main="Transformed GDP")
b.
fit a suitable ARIMA model to the transformed data using ARIMA();
US_economy_ts_transformed <- ts(US_economy$GDP_transformed, start=min(US_economy$Year), frequency=1)
fit_transformed <- auto.arima(US_economy_ts_transformed)
summary(fit_transformed)## Series: US_economy_ts_transformed
## ARIMA(1,1,0) with drift
##
## Coefficients:
## ar1 drift
## 0.4586 218.4542
## s.e. 0.1198 17.5690
##
## sigma^2 = 5488: log likelihood = -325.37
## AIC=656.74 AICc=657.19 BIC=662.87
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 1.540306 72.14009 54.34684 0.0007634341 0.4284988 0.2423587
## ACF1
## Training set -0.02316701checkresiduals(fit_transformed)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,0) with drift
## Q* = 3.8137, df = 9, p-value = 0.9232
##
## Model df: 1. Total lags used: 10c.
try some other plausible models by experimenting with the orders chosen;
Manually adjust the orders of the ARIMA model.
- Lower AIC/BIC values generally indicate a better model fit, as long as the residuals still look like white noise.
# ARIMA(2,1,0) model
fit_arima_210 <- Arima(US_economy_ts_transformed, order=c(2,1,0), include.drift=TRUE)
summary(fit_arima_210)## Series: US_economy_ts_transformed
## ARIMA(2,1,0) with drift
##
## Coefficients:
## ar1 ar2 drift
## 0.4554 0.0071 218.3917
## s.e. 0.1341 0.1352 17.7305
##
## sigma^2 = 5589: log likelihood = -325.37
## AIC=658.74 AICc=659.51 BIC=666.91
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 1.566054 72.13822 54.29574 0.001150807 0.427942 0.2421308
## ACF1
## Training set -0.02039424# ARIMA(1,1,2) model
fit_arima_112 <- Arima(US_economy_ts_transformed, order=c(1,1,2), include.drift=TRUE)
summary(fit_arima_112)## Series: US_economy_ts_transformed
## ARIMA(1,1,2) with drift
##
## Coefficients:
## ar1 ma1 ma2 drift
## 0.8296 -0.3892 -0.1937 214.7297
## s.e. 0.2039 0.2480 0.1586 23.5060
##
## sigma^2 = 5639: log likelihood = -325.11
## AIC=660.23 AICc=661.4 BIC=670.44
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 2.931548 71.78389 53.85649 0.01900945 0.4211449 0.240172
## ACF1
## Training set -0.01635229fit_arima_112_nodrift <- Arima(US_economy_ts_transformed, order=c(1,1,2), include.drift=FALSE)
summary(fit_arima_112_nodrift)## Series: US_economy_ts_transformed
## ARIMA(1,1,2)
##
## Coefficients:
## ar1 ma1 ma2
## 0.9944 -0.4986 -0.2366
## s.e. 0.0082 0.1314 0.1261
##
## sigma^2 = 5831: log likelihood = -327.73
## AIC=663.46 AICc=664.23 BIC=671.63
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 11.33192 73.67833 56.74494 0.1363926 0.4464121 0.253053
## ACF1
## Training set 0.006911175# ARIMA(1,1,0) model
summary(fit_transformed)## Series: US_economy_ts_transformed
## ARIMA(1,1,0) with drift
##
## Coefficients:
## ar1 drift
## 0.4586 218.4542
## s.e. 0.1198 17.5690
##
## sigma^2 = 5488: log likelihood = -325.37
## AIC=656.74 AICc=657.19 BIC=662.87
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 1.540306 72.14009 54.34684 0.0007634341 0.4284988 0.2423587
## ACF1
## Training set -0.02316701checkresiduals(fit_arima_210)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,1,0) with drift
## Q* = 3.838, df = 8, p-value = 0.8714
##
## Model df: 2. Total lags used: 10checkresiduals(fit_arima_112)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,2) with drift
## Q* = 3.7916, df = 7, p-value = 0.8034
##
## Model df: 3. Total lags used: 10checkresiduals(fit_arima_112_nodrift)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,2)
## Q* = 4.028, df = 7, p-value = 0.7766
##
## Model df: 3. Total lags used: 10checkresiduals(fit_transformed)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,0) with drift
## Q* = 3.8137, df = 9, p-value = 0.9232
##
## Model df: 1. Total lags used: 10d.
choose what you think is the best model and check the residual diagnostics;
Based on AIC and BIC, the ARIMA(1,1,0) model with drift is the best because it has the lowest AIC and BIC. The Root Mean Square Error (RMSE) and Mean Absolute Error (MAE) are also low compared to the other models, indicating a good fit to the data.
No trends were found in the residuals. Most of the points lie along the 45-degree line, which indicates that the residuals are normally distributed. Therefore, it is a relatively good model.
checkresiduals(fit_transformed)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,0) with drift
## Q* = 3.8137, df = 9, p-value = 0.9232
##
## Model df: 1. Total lags used: 10qqnorm(residuals(fit_transformed))
qqline(residuals(fit_transformed))
e.
produce forecasts of your fitted model. Do the forecasts look reasonable?
The forecast plot shows past trends continuing into the future, which appears to be a reasonable extrapolation of the time series data presented.
forecast_transformed <- forecast(fit_transformed, h=10)
plot(forecast_transformed, main="Forecast from ARIMA Model on Transformed Data")
f.
compare the results with what you would obtain using ETS() (with no transformation).
ARIMA’s AIC is lower than ETS’s and the ARIMA model contains a drift component that suggests it is explaining trends in the data. Additionally, because ETS has not been transformed, it may affect the results. Therefore, if we must compare, I think ARIMA is more suitable for the model.
ets_fit <- ets(US_economy$GDP)
summary(ets_fit)## ETS(M,A,N)
##
## Call:
## ets(y = US_economy$GDP)
##
## Smoothing parameters:
## alpha = 0.9991
## beta = 0.5012
##
## Initial states:
## l = 448093333333.703
## b = 64917355686.8708
##
## sigma: 0.026
##
## AIC AICc BIC
## 3190.787 3191.941 3201.089
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 20997049199 166906260796 104653662551 0.4843392 1.914575 0.3067446
## ACF1
## Training set 0.1532127ets_forecast <- forecast(ets_fit, h=10)
plot(ets_forecast, main="Forecast from ETS Model on Original US GDP Data")
plot(forecast_transformed, main="Forecast from ARIMA Model on Transformed Data")