786GroupProject.knit

Exploratory data analysis:

Write a brief description of the data set you have been given.
Visualise your time series (time plot, subseries plot, decomposition plot, etc) and discuss the features of the original time series.

training_data <- read.csv("qgdp_training.csv") %>% 
  select(c(Date, Retail.Trade)) %>% 
  mutate(Date = yearquarter(Date)) %>%
  as_tsibble(index = Date)
training_data %>% 
  autoplot(Retail.Trade) +
  labs(x = "Time",
       y = "Retail.Trade",
       title = "Time Plot")

training_data %>%
  gg_subseries(Retail.Trade) +
  labs(x = "Year",
       y = "Retail.Trade",
       title = "Subseries Plot")

training_data %>%
  model(stl = STL(Retail.Trade, robust = TRUE)) %>% 
  components() %>%
  autoplot()

2. ETS models: - Explore different ETS models for your data set. - Describe the methodology used to create a shortlist of appropriate candidate ETS models. You are expected to fit these both manually and automatically. - Select the ETS model that you believe has the best predictive ability. Explain your reasons for selecting this model. - Write down the fitted model equations.

# Fit the inferred ETS model
ets_manual_AMA <- training_data %>%
  model(ETS(Retail.Trade ~ error("A") + trend("M") + season("A")))
ets_manual_AAM <- training_data %>%
  model(ETS(Retail.Trade ~ error("A") + trend("A") + season("M")))
ets_manual_AAA <- training_data %>%
  model(ETS(Retail.Trade ~ error("A") + trend("A") + season("A")))
ets_manual_MAN <- training_data %>%
  model(ETS(Retail.Trade ~ error("M") + trend("A") + season("N")))
ets_manual_MAM <- training_data %>%
  model(ETS(Retail.Trade ~ error("M") + trend("A") + season("M")))
ets_manual_MMM <- training_data %>%
  model(ETS(Retail.Trade ~ error("M") + trend("M") + season("M")))
ets_auto <- training_data %>%
  model(ETS(Retail.Trade))

# Compare the AIC values of the models
aic_values <- bind_rows(
  ets_manual_AMA %>% glance() %>% mutate(Model = "ETS(AMA)"),
  ets_manual_AAM %>% glance() %>% mutate(Model = "ETS(AAM)"),
  ets_manual_AAA %>% glance() %>% mutate(Model = "ETS(AAA)"),
  ets_manual_MAN %>% glance() %>% mutate(Model = "ETS(MAN)"),
  ets_manual_MAM %>% glance() %>% mutate(Model = "ETS(MAM)"),
  ets_manual_MMM %>% glance() %>% mutate(Model = "ETS(MMM)"),
  ets_auto %>% glance() %>% mutate(Model = "ETS(Auto)")
) %>% arrange(AIC)

print(aic_values)

## # A tibble: 7 × 10
##   .model            sigma2 log_lik   AIC  AICc   BIC    MSE   AMSE     MAE Model
##   <chr>              <dbl>   <dbl> <dbl> <dbl> <dbl>  <dbl>  <dbl>   <dbl> <chr>
## 1 "ETS(Retail.Tra… 8.79e-4   -894. 1805. 1807. 1832.  5574.  6199.  0.0195 ETS(…
## 2 "ETS(Retail.Tra… 8.79e-4   -894. 1805. 1807. 1832.  5574.  6199.  0.0195 ETS(…
## 3 "ETS(Retail.Tra… 8.79e-4   -894. 1806. 1807. 1832.  5599.  6181.  0.0193 ETS(…
## 4 "ETS(Retail.Tra… 5.23e+3   -934. 1886. 1887. 1912.  4927.  5841. 41.4    ETS(…
## 5 "ETS(Retail.Tra… 6.73e+3   -951. 1921. 1922. 1947.  6338.  6966. 48.1    ETS(…
## 6 "ETS(Retail.Tra… 6.76e+3   -952. 1921. 1923. 1948.  6367.  6947. 49.1    ETS(…
## 7 "ETS(Retail.Tra… 6.52e-3  -1035. 2080. 2081. 2095. 29813. 28719.  0.0637 ETS(…

# Select the model with the lowest AIC value
best_model_name <- aic_values %>% slice(1) %>% pull(Model)
best_model <- switch(best_model_name,
  "ETS(AMA)" = ets_manual_AMA,
  "ETS(AAM)" = ets_manual_AAM,
  "ETS(AAA)" = ets_manual_AAA,
  "ETS(MAN)" = ets_manual_MAN,
  "ETS(MAM)" = ets_manual_MAM,
  "ETS(MMM)" = ets_manual_MMM,
  "ETS(Auto)" = ets_auto
)

# Print the best model equation
best_model_equation <- best_model %>%
  report()

## Series: Retail.Trade 
## Model: ETS(M,A,M) 
##   Smoothing parameters:
##     alpha = 0.403321 
##     beta  = 0.07330283 
##     gamma = 0.1493439 
## 
##   Initial states:
##      l[0]      b[0]      s[0]   s[-1]     s[-2]     s[-3]
##  1168.818 -2.777311 0.9594427 1.08881 0.9763337 0.9754133
## 
##   sigma^2:  9e-04
## 
##      AIC     AICc      BIC 
## 1805.476 1806.872 1831.887

print(best_model_equation)

## # A mable: 1 x 1
##   `ETS(Retail.Trade ~ error("M") + trend("A") + season("M"))`
##                                                       <model>
## 1                                                <ETS(M,A,M)>

ARIMA models:

Explore different ARIMA models for your data set.
Provide an explanation of any transformations and differencing required.
Describe the methodology used to create a shortlist of appropriate candidate ARIMA models. You are expected to fit these both manually and automatically.
Select the ARIMA model that you believe has the best predictive ability. Explain your reasons for selecting this model.
Write down the fitted model equation in backshift notation.

ARIMA_data <- training_data
ARIMA_data %>% autoplot(Retail.Trade) + labs(title = "ARIMA Retail Trade Forecasts", x = "Time", y = "Retail.Trade")

# Provide an explanation of log transformations and differencing required.
ARIMA_data %>% autoplot(log(Retail.Trade)) + labs(title = "ARIMA Retail Trade Forecasts", x = "Time", y = "Retail.Trade")

# These functions suggest we should apply only one seasonal difference(no first difference).
ARIMA_data %>%
  mutate(log_Retail_Trade = log(Retail.Trade)) %>% features(log_Retail_Trade, unitroot_nsdiffs)

## # A tibble: 1 × 1
##   nsdiffs
##     <int>
## 1       1

ARIMA_data %>% 
  mutate(d_log_Retail_Trade = difference(log(Retail.Trade), 12)) %>%
  features(d_log_Retail_Trade, unitroot_ndiffs)

## # A tibble: 1 × 1
##   ndiffs
##    <int>
## 1      0

# Check whether it is a stationary time series
ARIMA_data %>% autoplot(log(Retail.Trade) %>% difference(4))

## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_line()`).

# Plot the ACF/PACF of the differenced data and try to determine possible candidate models.
ARIMA_data %>% gg_tsdisplay(log(Retail.Trade)%>% difference(4),plot_type = "partial", lag_max = 16)

## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_line()`).

## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_point()`).

# Fit these both manually and automatically.
# Features of an AR(p) model the ACF is exponentially decaying or sinusoidal. There is a significant spike at lag p in the PACF, but none beyond lag p.
# Features of an MA(q) model there is a significant spike at lag q in the ACF, but none beyond lag q. The PACF is exponentially decaying or sinusoidal.
ARIMA_fit <- ARIMA_data %>%
  model(auto = ARIMA(log(Retail.Trade), stepwise = FALSE, approximation = FALSE),
        arima200210 = ARIMA(log(Retail.Trade) ~ pdq(2,0,0) + PDQ(2,1,0)),
        arima201210 = ARIMA(log(Retail.Trade) ~ pdq(2,0,1) + PDQ(2,1,0)),
        arima202210 = ARIMA(log(Retail.Trade) ~ pdq(2,0,2) + PDQ(2,1,0)),
        arima003210 = ARIMA(log(Retail.Trade) ~ pdq(0,0,3) + PDQ(2,1,0)))

# We select the best models by minimizing the AIC, AICc or BIC. We would prefer to use AICc.
glance(ARIMA_fit) %>% arrange(AICc) %>% select(.model:BIC)

## # A tibble: 5 × 6
##   .model        sigma2 log_lik   AIC  AICc   BIC
##   <chr>          <dbl>   <dbl> <dbl> <dbl> <dbl>
## 1 auto        0.000797    291. -566. -565. -543.
## 2 arima202210 0.000797    291. -566. -565. -543.
## 3 arima200210 0.000868    285. -558. -558. -541.
## 4 arima201210 0.000874    285. -557. -556. -536.
## 5 arima003210 0.000963    279. -544. -543. -524.

# Details of the optimal model are given
best_ARIMA_model <- ARIMA_fit %>% select(auto)
report(best_ARIMA_model)

## Series: Retail.Trade 
## Model: ARIMA(2,0,2)(2,1,0)[4] w/ drift 
## Transformation: log(Retail.Trade) 
## 
## Coefficients:
##          ar1      ar2      ma1     ma2     sar1     sar2  constant
##       1.8112  -0.8481  -1.4182  0.6660  -1.0612  -0.5813    0.0033
## s.e.  0.0790   0.0768   0.0878  0.0825   0.1067   0.1092    0.0006
## 
## sigma^2 estimated as 0.0007973:  log likelihood=291.24
## AIC=-566.49   AICc=-565.34   BIC=-543.24

Neural network autoregression (NNAR) models:

In Chapter 12 of the fpp3 online textbook, there is an advanced topic on the Neural Network Autogregressive (NNAR) model.
Read this chapter and any other materials that will help you understand the NNAR model.
In your own words, explain the NNAR model as if you were teaching it to your fellow classmates.
Fit an automatic NNAR model to your data set and report the fitted model.

fit_NNAR <- training_data %>%
  model(NNETAR(Retail.Trade)) %>%
  report()

## Series: Retail.Trade 
## Model: NNAR(7,1,4)[4] 
## 
## Average of 20 networks, each of which is
## a 7-4-1 network with 37 weights
## options were - linear output units 
## 
## sigma^2 estimated as 1506

Assumption checking:

Produce residual diagnostic plots for your three chosen models (ETS, ARIMA, NNAR).
Do not perform a Box-Pierce or Ljung-Box test.
Research the topic of surrogate data testing in time series, specifically, the Random Shuffle method. Explain this test for independence as if you were teaching it to your fellow classmates.
Test the independence of your residuals (or innovation residuals) using the surrogate_test.R code provided. Explain how the R code works.

source("surrogate_test.R")
# Assumption checking
residuals_MAM <- residuals(ets_manual_MAM, lag_max= 8)
gg_tsresiduals(ets_manual_MAM)

set.seed(1)
n_simulations <- 1000

# Calculate the autocorrelation of the original residuals
original_acf <- residuals_MAM %>% ACF(var = .resid) %>% pull(acf) %>% .[2]

## Warning: The `...` argument of `PACF()` is deprecated as of feasts 0.2.2.
## ℹ ACF variables should be passed to the `y` argument. If multiple variables are
##   to be used, specify them using `vars(...)`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

surrogate_acfs <- replicate(n_simulations, {
  shuffled_residuals <- sample(residuals_MAM$.resid)
  acf(shuffled_residuals, plot = FALSE)$acf[2]
})

# Draw a comparison chart of the autocorrelations of the original residuals and the substituted residuals
surrogate_acfs_df <- tibble(acf = surrogate_acfs)
p4 <- surrogate_acfs_df %>%
  ggplot(aes(x = acf)) +
  geom_histogram(bins = 30, fill = "blue", alpha = 0.7) +
  geom_vline(xintercept = original_acf, color = "red", linetype = "dashed") +
  ggtitle("Surrogate ACF Distribution vs Original ACF") +
  labs(x = "ACF", y = "Frequency") +
  annotate("text", x = original_acf, y = max(table(cut(surrogate_acfs, breaks = 30))), 
           label = paste("Original ACF =", round(original_acf, 3)), 
           vjust = 2, hjust = -0.2, color = "red")
print(p4)

# p-value
p_value <- mean(abs(surrogate_acfs) >= abs(original_acf))
cat("P-value for the independence test:", p_value, "\n")

## P-value for the independence test: 0.229

if (p_value < 0.05) {
  cat("The test rejects the null hypothesis of independence (p-value < 0.05).\n")
} else {
  cat("The test does not reject the null hypothesis of independence (p-value >= 0.05).\n")
}

## The test does not reject the null hypothesis of independence (p-value >= 0.05).

# test
lag <- 8
surrogate_result <- surrogate.test(residuals_MAM, lag = lag, N = 1000, test.stat = "ljung-box")
print(surrogate_result)

## $Q.null
##    [1] 18.0962283  6.5351514  8.0035683  7.3630809 15.9455589  5.7240429
##    [7]  8.6324950  6.3621187  7.4028805  8.3679081  7.5606072 10.1787658
##   [13]  5.3377269  4.3020013  5.8625599  9.2443101  7.9431016 13.7795924
##   [19] 13.3060276  6.4990901 19.1767417 11.6659182  8.3628524  7.1890070
##   [25]  6.8436774 11.6838662  6.6165318  6.2764710  9.8561999 16.5817268
##   [31]  8.0593623 16.6472260 15.0851826  8.4541531 19.5021623  2.2131561
##   [37]  3.9784133  3.6960737  5.6645615  9.9270868  6.7750724  4.3113328
##   [43]  7.2279088 15.2838349  6.8616059 10.7621216  5.9802215  8.3629710
##   [49]  2.2360168  7.7112856 11.1161213 11.1538827  4.9339150  8.0388893
##   [55]  5.3857367 18.9176767 10.8148693  8.2239087  6.0868985  3.4751827
##   [61] 10.1563459  2.7632292  6.1824930  5.5004664  6.8191912  0.7799284
##   [67]  5.0035683  2.5321771  7.4582687  6.9866877  2.0456570  9.4210976
##   [73]  7.3027381 12.9525481 11.3785247  7.8213380 12.5028158  3.4044164
##   [79]  1.7021423  3.3500496 11.2178435 15.3860473 10.1668701  5.8840210
##   [85]  4.8365290  6.4045939  3.2395502 11.7758211  8.9298911  5.2545349
##   [91]  9.6167855  4.1245964  6.6851024  4.1044025  5.3413993 14.8271630
##   [97]  8.0530142  3.6452892  5.3939405  4.7799950  4.0669119 10.5424532
##  [103]  4.9130061  3.5692074  6.0254287 13.6561823  3.7331883  2.3034363
##  [109]  4.4420636  8.2826288 15.6626277  6.8056216  5.8515644  9.9683176
##  [115]  1.7432988  5.3911115  7.8938870  8.5154372  5.5610292 16.2558877
##  [121] 12.9441523  6.4747321  5.9239898  4.4177577 15.2936138  6.7211514
##  [127]  5.4265638  6.0869824  5.0920697  9.4601696 14.6883401  6.0544757
##  [133] 12.7934501  6.4288159  4.5316716  7.5138607  8.0875490  7.6615067
##  [139]  9.5699951  7.6280434 10.6439385  3.7379548  6.9222815  3.3886577
##  [145] 17.0089162 13.2234329  9.2216033  3.9003982  4.3562047  5.8532793
##  [151]  9.3054952 12.5433354  2.1367211  6.9235586  5.9901355  2.5692264
##  [157]  6.3243190 11.2269005  7.7181619  7.6977147  5.1431349  5.3441131
##  [163]  9.7179410  5.9853483  7.1471340 10.3589336  4.5957507  1.9921200
##  [169]  2.4796151  7.1348627 17.1280721 11.5154565  6.3719811  7.2469577
##  [175]  2.4575170  9.6440305 10.7590556  6.1880472  4.0580805 15.7197964
##  [181] 17.8593842  4.6106203  3.3714319  8.9994372  4.5786594  4.1532291
##  [187]  5.5600615 12.3585372  3.8078389  6.9664544  4.9717200 13.9990740
##  [193]  4.2344088  6.3055088  6.1132193  8.0689465  7.5946950  6.9383572
##  [199]  2.2514233 11.2191946 12.9026244  9.4384553  7.4643169 12.8618057
##  [205]  7.8833642  8.1580272  4.9061982  5.2407787  5.4952156  1.4931072
##  [211]  3.7441339 14.0512524  4.0759511 15.1189953  3.5156546  5.6005361
##  [217]  4.2758273  7.2665302 10.2012771  3.0207694  8.4118971  9.8569523
##  [223]  2.5305213  6.6802740  2.9986631  3.5258314  3.0235915  6.4158322
##  [229]  8.0926522  2.5940641  2.6832044  4.8196431 15.1684033 11.5574380
##  [235]  8.4916643  5.1899112  7.4035277  9.0828285  9.2458602 12.0338820
##  [241]  5.6800608  3.9138268  6.3109395 16.6269610 14.9654139  4.9874374
##  [247]  7.6373315  8.2070167  6.9854298  7.2715958  8.2914360  5.1782278
##  [253]  9.1409660  7.7410888  8.0722868  4.1659676  6.8135773 12.0512063
##  [259]  5.4247103  2.1048139  4.9242748 11.6966146  8.5781741 13.8156123
##  [265]  3.7406800  4.8804096 10.6909687  6.6411673  9.1882765  2.2583368
##  [271]  9.9460676  5.3791705  2.7155477  3.8317861  6.6794559  2.7624142
##  [277]  9.4737485 14.9954886  9.1306750  5.3131008 10.5027262  5.2657532
##  [283]  2.4130029  5.3095869 11.0821960  9.6700521 12.5488944  5.2642345
##  [289]  8.0127975  2.2049697  6.3767667 13.9461292  3.6101945  5.4797327
##  [295]  7.8050753  4.8358548  2.2791794  6.5325165  6.9590740 10.5184774
##  [301]  3.6432354  6.3365125  7.2328058 16.7019821  7.5542608  8.2167181
##  [307]  2.7686207 12.0175592 13.4384232  5.4520975  3.8183690  2.4210633
##  [313]  6.1840443 18.7936220  5.1504930  1.1626337  3.5798035  5.5603631
##  [319]  6.2916209  5.4578928  3.7911097  5.7354374  8.4218915 13.9418459
##  [325] 13.5173893  5.8129213 14.4942616 33.2051123  3.7646962  4.7638019
##  [331]  2.2746721 12.8502013 12.4604303  3.1818538  6.1239340  5.4830942
##  [337] 13.3388420  4.7102046  6.1834709  6.3359004  5.7144638  9.1783700
##  [343] 14.7998354  3.4068915  3.6038636 12.5278012  3.0280976  6.8456799
##  [349]  8.5102319 19.3018105  3.4996842 11.7302389  5.8574957  5.2718183
##  [355]  6.8481731  9.1677720  6.4106725 17.7473657  4.7832267  6.2646024
##  [361]  2.8107920 17.7348671  8.9828220  5.7758095 17.3401181  6.6614741
##  [367]  7.6186149 29.1371519 26.4050195  4.6965593  5.0602013  4.4125514
##  [373]  6.5803666  3.9105832  5.6528247  9.4696352  4.1774715  4.8631002
##  [379]  7.6130452  5.9420309  6.2631171 15.6617800 16.8204720  5.9537662
##  [385]  5.7383184  8.3570963 11.2913825  6.4447292  9.1979124  3.7900889
##  [391]  5.0403165 11.7209055 16.0757733  2.7391030  3.2143386  8.8932375
##  [397]  7.7102500  3.3610031  6.7798475 12.6906488 10.5002532  5.4068371
##  [403]  6.6895366 19.4100489  5.7149070  2.2393808  9.9965598  7.6182455
##  [409] 11.4602807  4.7958130  6.1909739  4.4132658 10.4180705  1.5348362
##  [415]  6.4408544  5.8407405  9.6049315  7.3704055  5.6545593  7.9011807
##  [421]  7.3128098  8.6469385 12.2092460  8.6103537  8.4854668  6.0171974
##  [427] 14.1983263 10.5525809  6.8704418 11.4260127 13.0505260  8.4476638
##  [433]  4.0670475  6.8419564  3.1197157  6.9256068  8.5124237 12.6024808
##  [439]  6.6007312  3.0757189  6.1312015  2.3977057  9.7673023  9.7266728
##  [445] 18.1612770  2.9102357  4.4730823  6.3494901  6.6387122 10.7245882
##  [451]  6.5888441  5.2477835  8.0040222  5.9513282  4.4167526  5.3019028
##  [457]  7.2859041  3.4442679  3.4792272  6.3264444  5.7417279 15.6155173
##  [463]  5.9465235 15.8971386 22.9633987  4.5557166 15.8419111 13.3068538
##  [469] 11.2003466  1.8798189  8.9193137  4.4276997  6.1078997 12.5638319
##  [475]  6.5083559 16.9004106 13.3778173 14.2788166  6.7633311  5.6608445
##  [481] 10.8363456 12.4561686  5.3984119  5.7958634  2.2955564  7.3988748
##  [487] 10.1839897  4.3288609  8.0314664  5.4387117  3.6840873  6.4934137
##  [493]  9.4747733  2.1329658  5.9161688  7.5040409  6.6126352  3.7725955
##  [499]  8.2614344 11.3666876  1.6040005 10.6455058  4.5350840  7.2580296
##  [505]  4.5608227 14.4242868 28.7190334  8.6015252 20.7468546  5.6728553
##  [511]  8.4311191  5.3760825  4.9162451 16.6701121  6.1022206  4.8093999
##  [517] 10.5481424  6.4909368  5.1443097  6.4153081  2.4973241 10.5342507
##  [523]  6.2364674  6.1631722  8.0168521  8.2859198  4.4729281 14.7842509
##  [529] 10.4442556  6.8282213  6.4583234  6.1230606  5.7591248  8.5003479
##  [535]  6.6184067  4.8668416  8.0603600  3.7676160  4.1759245  9.3991181
##  [541] 10.8663248 18.4000548  9.1698314 11.7295888  7.3141199  4.7006637
##  [547]  2.6068516  6.7590682  4.9598845  5.1650601  6.6450879  4.6493903
##  [553]  8.9666053  3.3178882  3.0043837  6.8551293  5.6125940  6.4970000
##  [559] 10.3916786 18.5712242  7.4738426  2.7066667 10.4268987  1.6646876
##  [565] 19.9945967 31.0177949  6.6564048  2.3622768 16.6562922  9.0370067
##  [571] 13.4669284 16.9413066  5.8059878  8.2814152 12.2539150  3.3554559
##  [577]  5.1655039  8.8298509  6.0815615 11.6723607  2.9496138  5.0237194
##  [583]  6.8218305  7.9169085 10.2541758 10.6445469  8.1376603  6.0384652
##  [589] 12.3410171  2.1329040  9.2341026 19.1934068 15.7237934 11.8612171
##  [595]  9.9798509  2.5870965  4.1964269  7.4260814 11.8234955  4.6542679
##  [601]  8.4657552 19.2518339  6.8279498  9.1129357  3.4432835  2.5488631
##  [607]  7.6114201  3.4570444  7.9590588  7.6382228  6.8846849 14.8238807
##  [613]  9.0252214  6.4623694  8.7091417  6.3387635 11.8160383  9.0179988
##  [619]  5.8002882  4.8599809  2.7056718  9.8823230  2.6092786 12.5415714
##  [625]  9.1016868 10.8545729  6.1715584  4.4251478  6.5783462  3.7771058
##  [631]  9.7965840  4.3630845  4.4037981 14.5347830  3.3796167  5.9021094
##  [637]  6.0431152  1.3979155  7.3305021  8.1115896  5.1848827  3.8165197
##  [643]  9.2934723 14.0662047  5.6352320 10.8946861 10.6655931 10.0138625
##  [649]  8.2380409  8.8505899 12.2253464  8.7223064  3.0970043  3.7074576
##  [655]  8.5467105 12.2057286  7.4510091  7.2473457 10.4969926  3.0208423
##  [661]  8.3907844  7.7218748  9.4018672  9.8146448  2.4111216 25.6963888
##  [667]  6.9604120  2.1189158  5.3704310  3.9715186  3.4595027  9.7189244
##  [673] 19.5347740  1.9699662  7.6284471  3.7391845  3.6372308  6.9738454
##  [679]  4.0986279  3.5777381  7.8971137  9.2269618  8.3629271  4.0366852
##  [685]  6.1784058  3.9495293  2.2623940 10.9017181  7.2100662  4.0527083
##  [691]  4.9120550 15.4335039 15.9113636  9.9627528  1.8649064  7.1027277
##  [697] 13.7161008  5.7382384  4.1176317 14.2254662  4.2693841  8.9101803
##  [703]  3.4970762  7.3060321 10.3342646  8.0159333 11.3532171  4.6988353
##  [709] 11.7178167  6.1168924  9.7237338  3.1093827  5.3467712 10.4403220
##  [715]  7.4518885 14.2995664  6.1482303  7.1156651  1.7931234  3.0382801
##  [721]  8.2411451  9.5069819  9.3495311  6.1788460 14.9391673 13.0536260
##  [727]  2.3561440  6.6007807  8.7128183  6.8978610  5.8877487  8.9793699
##  [733]  5.5115753  4.6127668 13.1106627  7.5208224  9.0257874  5.7852335
##  [739]  3.1755462 18.1230076  7.5106084 10.4441964  5.6764798 18.2530610
##  [745]  6.1868809  3.6643527  2.7518226  7.1961017  1.7440889 10.3041088
##  [751] 12.7161003  4.5111864 11.3709188 12.5091952  9.6995443 12.2770674
##  [757]  8.9013120 19.7351668  5.3411319 10.8837528  6.8374538  7.7884612
##  [763]  4.4513935  6.6269498  7.0474580  8.8504930  7.1789285  5.1224553
##  [769]  6.0591150 13.4765153  4.8776338  7.2884572  4.8219311  5.8527450
##  [775]  6.3134333  8.2640155  6.5988859  6.7347422  6.7080213  7.6583103
##  [781]  7.6931102  3.3939630  4.7450717  2.8436476  2.4400475  5.6520697
##  [787]  6.8850586  6.5030285  7.3035790  8.3631844  5.3844047  5.8324877
##  [793] 10.8675423  9.4206022  8.7141899  5.9967679  6.3065290  5.4444036
##  [799] 16.0192296  5.7884370 10.3551256  8.8454794  6.1462013  5.3277164
##  [805]  2.7580071  9.4198319  2.2442040  2.1641363 13.1673262  3.2328812
##  [811]  6.1930686  7.2174756 13.4742984  7.1152355 12.5190869  7.0641343
##  [817]  8.4842349  6.2021071 17.8327158  7.0250266  6.9584258  7.8816496
##  [823]  9.5188189  2.6760197  2.7495445  9.6204229  1.8624547  8.1289014
##  [829]  6.5762309  1.9031416  5.1346055  8.3031767 12.3252580  1.8483222
##  [835]  9.0940055  3.5942379  5.9426680  9.5263982  2.5625598  5.5249911
##  [841] 12.0202469  4.8922084 18.1842121  7.8952417  6.3054738  6.3721081
##  [847]  4.9670847  8.9887944 10.5244962  4.3215345 12.3312423  8.2849582
##  [853]  9.8105640  5.3495648 14.7476750  6.3596382  5.7504187 10.2299726
##  [859] 11.7434591  7.9525503 13.4432370 16.6393041  3.6987209  5.0486351
##  [865]  8.8306367  4.0989996  9.5125737  9.2456189  9.0533041  3.3460141
##  [871]  6.5244347  1.7799722 12.0096623 12.1643701  9.7412834  4.8991521
##  [877]  1.6509129  6.1138951  5.3298830 12.2423231  6.9917234 10.6039746
##  [883] 13.6626084 12.7139333 13.5027450  5.8304815  5.1326515 23.8219325
##  [889]  4.5135022  6.0181111  4.4990182  6.9761921  2.6824118  7.8470692
##  [895] 10.8941140  6.3230242  5.8031920  9.3139970 22.9958731  7.9642620
##  [901]  5.7888184  5.9317222  5.3593934  6.6180294  5.3890608 22.8200678
##  [907]  7.4521479 16.6929986  7.3979342  7.1035777  6.9214012  3.2939702
##  [913]  5.3525557  9.5871664 14.5885288  7.5352865  6.9516791 10.7272755
##  [919] 11.2234926  9.6038694  4.6044848  4.0011486  6.6972489  3.3458062
##  [925]  5.5121736  9.2815259  5.1733575  5.3990303  5.3363486  8.1384180
##  [931]  8.3303069  8.4076632  4.3935681  9.2166646 11.3740777  7.0261568
##  [937] 14.0398291  1.5791570  6.2258880  9.9492834  5.5016417  8.8067824
##  [943]  9.9886831  5.6016347 11.4465707  7.1387231  3.9977572  7.4078854
##  [949]  2.4033466  6.9361471 10.1489379  4.4303257  5.3844371 10.1747864
##  [955]  5.5154930 10.0991099  8.1851999 14.9820277  4.6622828 11.4979388
##  [961]  0.7178135  5.2117874 11.4419922 10.1190588  4.3959519  5.1078008
##  [967]  4.1136680  7.7896916  1.7443548  8.4006172  5.2508425  6.0991514
##  [973]  8.0779902  4.9816326  5.0704310  7.1488888 12.2911934  0.5524640
##  [979]  2.5230461  4.2740712 10.8096540  6.0332270  6.3591128  6.7062187
##  [985]  5.2250274  9.7362521  6.4530712  9.8022328 15.0873036  5.6418964
##  [991]  9.7157304 13.3784091  6.4606074  2.2929056  8.3332560  3.4081244
##  [997]  9.1303861  4.3904879 11.9555948  4.2902343
## 
## $Q.obs
## [1] 6.425214
## 
## $test.stat
## [1] "ljung-box"
## 
## $p.value
## [1] 0.562
## 
## attr(,"class")
## [1] "surrogate"

plot.surrogate(surrogate_result)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#Check the residuals from your chosen model by plotting the ACF of the residuals.
best_ARIMA_model %>% gg_tsresiduals()

#Research the topic of surrogate data testing in time series
ARIMA_residuals <- residuals(best_ARIMA_model)

# 设定随机种子和模拟次数
set.seed(123)
n_simulations <- 1000

# 计算原始残差的自相关性
original_acf <- ARIMA_residuals %>% ACF(var = .resid) %>% pull(acf) %>% .[2]
# 生成替代数据并计算自相关性
surrogate_acfs <- replicate(n_simulations, {
  shuffled_residuals <- sample(ARIMA_residuals$.resid)
  acf(shuffled_residuals, plot = FALSE)$acf[2]
})

# 绘制原始残差和替代残差自相关性的比较图
surrogate_acfs_df <- tibble(acf = surrogate_acfs)
p4 <- surrogate_acfs_df %>%
  ggplot(aes(x = acf)) +
  geom_histogram(bins = 30, fill = "blue", alpha = 0.7) +
  geom_vline(xintercept = original_acf, color = "red", linetype = "dashed") +
  ggtitle("Surrogate ACF Distribution vs Original ACF") +
  labs(x = "ACF", y = "Frequency") +
  annotate("text", x = original_acf, y = max(table(cut(surrogate_acfs, breaks = 30))), 
           label = paste("Original ACF =", round(original_acf, 3)), 
           vjust = 2, hjust = -0.2, color = "red")
print(p4)

# 执行替代数据测试
ARIMA_lag <- 4 # 季节性数据滞后数(4的倍数)，可以根据具体情况调整
ARIMA_surrogate_result <- surrogate.test(ARIMA_residuals, lag = ARIMA_lag, N = 1000, test.stat = "ljung-box")
# 打印测试结果
print(ARIMA_surrogate_result)

## $Q.null
##    [1]  4.3512470  4.6515483  1.0655355  1.5668542  1.9257315  0.9911097
##    [7]  3.7801453  1.5566794  1.8408436  8.0799906  1.6084747  5.8441950
##   [13]  7.5031396  1.9552758  3.7320795  5.3281852  3.1703779  0.5074368
##   [19]  0.5582775  6.9579191  4.3751798  2.3001191  7.7792025  3.5596533
##   [25]  0.8343854  1.1030057  2.7658822  7.3056049  4.0597829  1.5035809
##   [31]  3.8592685  5.9403171  4.2848471  1.9283741 18.0237784  6.1441758
##   [37]  4.1566997  0.8777840  3.2899139  3.0037303  2.4577697  1.6316644
##   [43]  4.1299822  2.3683170  2.4369513  1.8257033  1.5862973  1.3945866
##   [49]  1.7444905  4.1328435  2.9552959  4.6994469  1.0446321  3.8263765
##   [55]  1.3016768  8.5266077 16.1810417  2.2652424  1.8729550  5.3299138
##   [61]  4.9576743  0.4600090  3.6752283  1.1706295  2.0393798  0.5172364
##   [67]  3.5955227  3.5654287  0.9827742  2.7440890  4.6703229  2.5253501
##   [73]  4.1441286  3.1603444  4.1682864  1.7548068  6.2571740  7.4493797
##   [79]  0.9197354  2.4567766  3.5064743  1.4363439  4.7401112  2.9635328
##   [85]  4.3678792  1.2634755  1.7897546  3.2302091  1.0480815  5.0310718
##   [91]  3.5254901  3.1135304  2.4092738  5.8579760  2.2088950  2.6115096
##   [97]  1.5836308  1.3737414  2.6023317  2.8402128  2.8395054  2.6885288
##  [103]  5.3285979  5.6839574  3.0172409  1.8656984  1.8424168  0.4181691
##  [109]  7.5955784 11.0944258  3.8779633  1.7125850  5.6286622  1.4555889
##  [115]  2.3548432  7.1217044  1.8307872  3.5634394  2.3915680  3.2933053
##  [121]  3.8201452  2.1606987  5.2110324  2.2953863  9.2555713  2.3070393
##  [127]  1.1969597  2.9603241  1.1423842  2.2009670  0.8861013 12.5533946
##  [133]  8.7418395  1.0531494  2.3647175  4.8509034  2.2645560  5.6210855
##  [139]  5.2552803  2.9990077  4.5709928  2.3122491  1.8983647  4.2960898
##  [145]  3.6029933  9.3536010  6.2542973  6.5246430  2.8166923  3.7594339
##  [151] 12.3912295  6.7064978  1.1377195  3.0035743  6.2671722  2.3441035
##  [157]  0.8627513  1.7551206  2.7107960 10.3143786  3.7868175  8.5294869
##  [163]  2.1565962  3.0760769  2.6279220  9.0898148  4.0610727  1.7766051
##  [169]  1.1115734  2.7744130  3.4993563  8.7032556  1.3482642  0.5781956
##  [175]  2.0561427  1.7351096  3.3334407  0.8588720  2.0789556  4.0051841
##  [181]  4.4762833  3.4716667  5.0323994  6.3704264  0.2299906  1.1613053
##  [187]  2.0935745  2.6106620  4.3832066  5.2159028  6.3701841  1.6979777
##  [193]  3.3531645  2.6130538  0.6652803  2.0072781  4.1930202 10.7453099
##  [199]  2.4322248  2.2615068  4.3266925  0.9530071  4.7352319  5.8083846
##  [205]  3.9972606  4.3271808  2.0351578  3.6350732  2.9358384  2.7565714
##  [211] 15.0348094  3.0044445  0.8819681 10.6109782  1.6837102  4.5811104
##  [217]  0.3510692  6.3338106  2.7152621  7.1541864  3.1694578  1.0463878
##  [223]  2.3042477  0.9780226  3.0379344  1.3965158  7.5543629  6.4165711
##  [229]  1.7800957  1.4062827  2.0782373  4.6517690  4.6528383  2.4783392
##  [235]  4.5258792  5.3285193  2.5416945  6.8368838  3.1618900  0.9624985
##  [241]  2.0022200  0.4627082  5.7382130  1.0972147  5.7414501  1.3952010
##  [247]  3.1760717  2.2241498  2.3941935  2.6881361  7.2469580  2.6912268
##  [253]  5.6172367  2.5302130  0.6690999  3.1842190  1.9452214  4.2407283
##  [259]  5.4729861  1.2934069  1.0276023  1.2607830  3.9046021  4.0016654
##  [265]  3.3275768  1.7029073  3.6795666  1.9898587 10.0228083  1.5168122
##  [271]  3.8008326  1.9533498  0.6580758  1.2147418  4.7271241  2.4967357
##  [277]  1.9360167  5.4720458  1.3414815  4.3445204  2.7349547  7.5728116
##  [283]  6.8989964  2.0458554  4.3681628  2.6282361  3.5296729  3.4736375
##  [289]  4.5014716  5.0452533  2.2216000  4.1096966  1.7149170  7.6080148
##  [295]  2.8775707  2.2616309  5.4691511  0.7409734  6.6758210  6.8474273
##  [301]  3.7723515  1.5935461  6.0085697  3.6162368  6.4748207  2.9580605
##  [307]  4.2887184  3.1198454  4.7425091  3.2941266  1.3675534  6.3661267
##  [313] 11.9141512  2.1648666  1.5173213  5.3457504  6.8633895 14.3173141
##  [319]  3.1993925  4.6735929  2.5869892  0.5155373  0.6859438  4.3403966
##  [325]  1.9506973  3.6078415  0.6538979  0.4248347  0.7743627  0.9167813
##  [331]  1.8028062  2.2914101  0.9587937  5.2658906  4.6652836  1.6331295
##  [337]  5.7949234  3.0989878  3.0124466 10.5315988  6.0273103  0.3845795
##  [343]  0.7787747  2.5053150  7.5028008  5.9617479  3.9285642  0.7994756
##  [349]  1.8823980  6.7387413  0.8429960  2.1740117  2.9493610  6.8309227
##  [355]  2.4491105  4.0565812  6.5038447  0.7209558  5.0448572  2.3741524
##  [361]  8.3092791  1.6005182  3.8950206  1.1753386  6.3051482  1.1151605
##  [367]  2.1471044  2.3226359  6.1450647  0.8584846  1.6661702  6.5017982
##  [373]  2.1634338  0.9064526  8.3973592  2.3358894  0.7909780  3.5903343
##  [379]  2.5940717  6.5820003  3.6156274  3.3685173  9.4496100  4.5562490
##  [385]  3.4143611  9.7090796  7.8146384  2.9377766  1.8817487 13.6287685
##  [391]  2.3588614  0.9328620  0.9606466  2.0898132  2.1235918  3.2144536
##  [397]  3.3640956  1.3387145  2.8128692  1.7170486  4.5307577  1.9756577
##  [403]  2.0164880  6.0934951  3.7376532  2.3941253  0.7987405  2.8521905
##  [409]  0.8013506  3.2432819  3.4054934  0.5965406 14.1379257  1.8335138
##  [415]  5.5249300  2.0866516  3.0898027 10.5828590  1.4155889  8.2714101
##  [421]  6.1958340  7.0544347  2.5028079  9.2260467  0.4888381  2.7912147
##  [427]  6.4813133  2.3874168  4.2901422  2.6313083  1.7715202  1.3782760
##  [433]  1.1602679  1.1571964  1.8522128  4.7392081  0.6267250  1.6627604
##  [439]  0.2609937 11.2935354  1.8511243  3.5886016  4.0581467  2.2769822
##  [445]  4.4342454  1.5123222  4.1318746  6.1304285  7.4081815  3.2964124
##  [451]  3.0103266  1.4482970  3.5655515  2.3342891  0.4114660  1.1001863
##  [457]  4.5366662  0.7425048  2.0416590  9.9767864  1.1468667  0.6412213
##  [463]  5.4124629  6.2195955  0.4097709  3.3120082  3.5642020 15.9098007
##  [469]  2.6604959  3.1667389  8.2157232  1.6064730  0.9699043  1.3148764
##  [475]  2.1901195  1.0558856  1.7027656  1.4242450  2.4422107  2.3724941
##  [481]  1.2017359  5.0178190  1.6036254  0.7706131  0.8148412  9.2651943
##  [487]  3.3728851  4.8850348  1.3047074  6.0797458  8.4549323  1.7890948
##  [493]  5.0724942  2.6759616  3.7102514  3.5321904  0.7733077  3.5964909
##  [499]  2.8421759  3.5530140  1.5745044  2.3224638  3.5178084  0.7912169
##  [505]  5.0077211  7.6441193  0.6836246  1.7755444  0.4571456  1.0608227
##  [511]  5.5010702  3.7187863  1.1406706  5.8953855  2.0542798  3.6170504
##  [517]  3.8789151  2.7384281  3.4300923  3.8169229  3.0847708  2.2009217
##  [523]  1.8809969  1.1058624  3.9748512  0.4575317  1.1599557  3.0100356
##  [529]  5.3724209  0.9696959  7.7980254  3.0478152  2.7465115  1.0957131
##  [535]  2.2394885  0.6974822  6.4749928  3.6103594 10.3112331  2.4310025
##  [541]  3.2861574  8.2737066  4.1151440  7.6509965  4.3227642  0.4281874
##  [547]  7.9786214  3.1823655  1.1927440  5.2766521  1.6731111  3.7200996
##  [553] 10.7092195  2.7009229  2.2027342  4.0596098  4.5402685  6.6108494
##  [559]  4.2758057  1.3807129  5.9773263  1.4842257  4.2362801  0.5415392
##  [565]  1.0347592  2.6993904  1.3342247  5.6565546  2.6361771  1.7295338
##  [571]  8.7672402  0.3552467  1.6321302  6.2221559  2.0464021  1.5136556
##  [577]  7.1703607  5.1149027  9.0172413  2.5663340  5.0778093  4.2330652
##  [583]  4.5837956  0.9254712  6.0911838  0.7762740  0.1615619  1.1420263
##  [589]  1.0373328  4.9900429  3.5802404  2.7928524  6.6326869  4.0650609
##  [595]  5.7176577  4.8796693  2.0727903  3.1918143 11.7329119  6.6064809
##  [601]  5.7021396  1.0134695  2.1840422  2.9832051  0.4330369  0.6762794
##  [607]  4.0699301  2.4822447  0.2655137  2.1441362  5.2767525  2.1660966
##  [613]  6.8073871  3.3076201  0.4977814  9.4002207  3.0875105 12.4576568
##  [619]  9.9323673  9.6189543  9.0748344  1.2526203  1.9206385  2.4795314
##  [625]  3.4535819  8.1133732  2.8155741  1.0699267 12.7603546  4.4737278
##  [631]  1.9547459  1.5889596  1.7174263  4.7932591  0.3092089  1.1692323
##  [637]  1.9075096  3.9391457  0.4547548  3.4345976  2.6266328  0.5484630
##  [643]  2.1005171  5.2983055  3.7716367  6.0898370  0.9522291  4.9496578
##  [649] 10.9117504 19.3012265  4.1718104  5.1439887  4.4427157  4.0958288
##  [655]  1.2627836  2.3329072  3.3303099  3.6472077  6.5470079  8.6454379
##  [661]  1.3717042  2.8092176  4.8422609  0.7965382  1.8931114  2.7868593
##  [667]  1.8527430  1.0670834  3.1105726  5.3993864  1.1799715  7.2785339
##  [673]  8.2149527  0.5233001  0.8007711  2.5080000  2.9982723  3.2954703
##  [679]  0.4073371  4.7557433  1.3216929  5.4819341  3.7263243  6.5292183
##  [685]  2.5827166  4.9500634  2.1269749  1.4601674  1.7478226  1.3860200
##  [691]  3.3518197  4.4634185  5.7016648  2.1447633  3.2482388  3.6190504
##  [697]  1.3150526  2.2148857  2.5981729  2.1063941  3.0105919  5.0119917
##  [703]  1.8340315  6.0345080  2.0797510  2.7128132  1.4818158  2.9447783
##  [709]  5.0600490  6.4624976  4.7538225 10.8306040  1.9328208  8.2556712
##  [715]  3.4964947  3.0218785  4.9081487  8.5542131  4.0986871  2.4271350
##  [721]  8.1415015  1.6882982  1.0915535  4.7561410  2.7499763 10.9211755
##  [727]  6.8563922  2.0094868  3.3258816  1.3045173  0.4287479  0.8191792
##  [733]  7.0093012  1.8234046  1.0979995  2.0175425  9.2374859  4.1312201
##  [739]  3.5832633 10.2310023  3.1426169  4.8144730  0.3611725  7.8143818
##  [745]  0.6828754  5.4299603  4.5236668  1.2301733  1.1461665  7.6901014
##  [751]  4.5259455  2.7535195  0.7933454  1.3778988  0.4014984  1.4748347
##  [757]  3.3919374  2.4578835  1.0858211  4.3615118  3.0822499 22.4301682
##  [763]  5.1082168  6.8510822  0.6642824  0.1641762 12.4897802  6.2958855
##  [769]  1.5501425  8.2211713  1.9045148  2.0247375  2.6951515  2.8552280
##  [775] 14.2605437  0.5371295  0.7617078  0.7565330  2.4124368  3.2295618
##  [781]  1.6880185  2.4703687  1.1897908  1.9410128  1.6743481  2.9734682
##  [787]  3.4268269  3.2629787  1.7189768  9.9906391  1.5398992  2.5797073
##  [793]  6.7657740  0.9914949  5.8052516  3.5685258  4.0320559  1.8208692
##  [799]  5.1984595  2.8595768  1.0260485  0.7141244  0.5443696  5.0110820
##  [805]  2.0606248  3.9916474  4.5837851  0.3587166  2.4549843  3.6236251
##  [811]  2.2184421  2.9196431  8.2117417  4.0041260  2.6868871  6.9040042
##  [817]  1.9518156  6.5174472  3.3757002  1.4960856  0.3619813  0.3251998
##  [823] 10.2848029  5.6198798  1.5106910  1.9747147  6.6582158  0.4881114
##  [829]  2.7392539  4.0661796  7.8908890  2.8628330  4.0010579  2.0057380
##  [835]  3.6262124  0.7144778  1.8177542  2.2526019  1.0308983  9.9674704
##  [841]  5.7100363  6.0646799  1.4276016 18.3334719  1.4880171  1.9691710
##  [847]  3.1066186  0.3011860  7.6095695  7.7182955  1.1603371  3.4362245
##  [853]  0.1758972 12.0040202  1.4730332  2.1727069  3.8162309  4.3571630
##  [859] 11.1635463  7.9485044  3.9691983  0.6482233  2.3306552  1.1117067
##  [865]  3.1321219  4.5746740  4.3890213  8.3771037  4.1858666  2.2480443
##  [871]  9.3708746  2.1097983  1.4228970  1.5509939  2.8989285  3.0958830
##  [877]  1.4184168  1.6376439  5.3037760  2.2123104 10.2622100  2.0267441
##  [883]  8.3732166  4.4557241  1.7510382  4.7413522  3.4892172  4.8856291
##  [889]  1.9373760  9.3693503  3.2604745  3.8032615  2.7921334  4.3416330
##  [895]  2.7580300  4.3436083  0.8913323  4.9380691  2.9193647  2.1753892
##  [901] 18.5266096  2.1447400  6.3270111  6.7433129  1.6474337  0.4672436
##  [907]  2.1427868  2.6913915  4.3938240  1.5314460  3.2031420  6.8878737
##  [913]  5.1241813  1.7694882  2.5114403  4.8894669  3.1490704  6.0424641
##  [919]  1.8178409  2.2013396 14.3109999  1.4688074  4.6937865  5.3193870
##  [925]  7.8372586  3.6182082  3.0474072  1.6239637  4.9158861  8.1754045
##  [931]  2.1491574  1.0913464  1.6585622  4.5982018  3.0823803  3.1232266
##  [937]  5.7243608  1.6984555 14.4465808  9.8792177  2.6199907  4.2030894
##  [943]  8.3523336  7.5556558  9.6435263  3.4528848  3.5603256  2.5377215
##  [949] 10.3188466  2.2118069  6.1932552  2.4915332  1.5822389  7.8223285
##  [955]  1.7909272  4.9408365  4.0314167  0.1401512  0.9409727  1.1208434
##  [961]  2.8952451  8.0988450  1.5518243  6.8098695  2.8745852  5.4334453
##  [967]  5.0219864  5.7670795  3.0983250  0.2810242  4.5299071  1.2601695
##  [973]  5.2954168 10.1707500  4.5939303 11.7543588  7.7282508  5.6723105
##  [979]  2.3343903  9.1224940  8.0086740  0.7332820  5.6497652  5.4747428
##  [985]  5.2997380  3.0622495  5.0043301  1.3861727  1.3379506  3.6699265
##  [991]  1.0281607  2.9061707  4.5975379  2.0618171 11.9210816  2.3489540
##  [997]  3.5095319  4.4525004  3.7171670  6.5265265
## 
## $Q.obs
## [1] 0.9412394
## 
## $test.stat
## [1] "ljung-box"
## 
## $p.value
## [1] 0.903
## 
## attr(,"class")
## [1] "surrogate"

# 绘制替代数据检验统计量分布图
plot.surrogate(ARIMA_surrogate_result)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# Residual diagnostic plots
fit_NNAR %>% gg_tsresiduals()

## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_line()`).

## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: Removed 7 rows containing non-finite outside the scale range
## (`stat_bin()`).

# Test the independence of your residuals
s = surrogate.test(na.omit(residuals(fit_NNAR)), lag=8, N = 1000, test.stat = "ljung-box")
s$p.value

## [1] 0.015

### The p-value > 0.05, we do not reject residuals from the model resemble a white noise.
s %>% 
  plot.surrogate() + 
  theme_bw() +
  labs(title = "Surrogate data test with Box-Pierce statistic")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Forecasting:

For your three chosen models, produce forecasts (point-forecasts and prediction intervals) for the h=8 quarters up to and including 2023 Q4.
Explain how the prediction intervals are calculated for your chosen models.
When the test set becomes available on Canvas 48-hours before the due time, compare your forecasts with the test set by evaluating model accuracy.
Write a small discussion on how well your three models performed at forecasting.

# Forecasting
h <- 8  
forecast_MAM <- ets_manual_MAM %>%
  forecast(h = h)

print(forecast_MAM)

## # A fable: 8 x 4 [1Q]
## # Key:     .model [1]
##   .model                                               Date   Retail.Trade .mean
##   <chr>                                               <qtr>         <dist> <dbl>
## 1 "ETS(Retail.Trade ~ error(\"M\") + trend(\"A\") … 2022 Q1 N(3568, 11188) 3568.
## 2 "ETS(Retail.Trade ~ error(\"M\") + trend(\"A\") … 2022 Q2 N(3555, 13565) 3555.
## 3 "ETS(Retail.Trade ~ error(\"M\") + trend(\"A\") … 2022 Q3 N(3598, 17157) 3598.
## 4 "ETS(Retail.Trade ~ error(\"M\") + trend(\"A\") … 2022 Q4 N(4286, 30148) 4286.
## 5 "ETS(Retail.Trade ~ error(\"M\") + trend(\"A\") … 2023 Q1 N(3761, 31402) 3761.
## 6 "ETS(Retail.Trade ~ error(\"M\") + trend(\"A\") … 2023 Q2 N(3745, 37568) 3745.
## 7 "ETS(Retail.Trade ~ error(\"M\") + trend(\"A\") … 2023 Q3 N(3788, 46129) 3788.
## 8 "ETS(Retail.Trade ~ error(\"M\") + trend(\"A\") … 2023 Q4 N(4509, 77985) 4509.

full_data <- read.csv("qgdp_full.csv") %>% 
  select(c(Date, Retail.Trade)) %>% 
  mutate(Date = yearquarter(Date)) %>%
  as_tsibble(index = Date)

# plot
autoplot(full_data, Retail.Trade) +
  autolayer(forecast_MAM, level = NULL) +
  labs(title = "Retail Trade Forecasts", x = "Time", y = "Retail.Trade") +
  guides(colour = guide_legend(title = "Models"))

autoplot(full_data, Retail.Trade) +
  autolayer(forecast_MAM) +
  labs(title = "Retail Trade Forecasts", x = "Time", y = "Retail.Trade") +
  guides(colour = guide_legend(title = "Models"))

accuracy(forecast_MAM, full_data) %>% 
  select(.model, RMSE, MAE, MAPE, MASE)

## # A tibble: 1 × 5
##   .model                                                  RMSE   MAE  MAPE  MASE
##   <chr>                                                  <dbl> <dbl> <dbl> <dbl>
## 1 "ETS(Retail.Trade ~ error(\"M\") + trend(\"A\") + sea…  457.  393.  11.3  4.27

# produce forecasts (point-forecasts and prediction intervals) for the h=8 quarters up to and including 2023 Q4.
best_ARIMA_model %>% 
  forecast(h = 8) %>%
  autoplot(ARIMA_data, level = NULL) + 
  labs(title = "ARIMA Retail Trade Forecasts", x = "Time", y = "Retail.Trade")

best_ARIMA_model %>% 
  forecast(h = 8) %>%
  autoplot(ARIMA_data) + 
  labs(title = "ARIMA Retail Trade Forecasts", x = "Time", y = "Retail.Trade")

ARIMA_forecast_data <- forecast(best_ARIMA_model, h = 8)

ARIMA_actual_values <- full_data %>% 
  tail(8) %>% 
  pull(Retail.Trade)
ARIMA_predicted_values <- ARIMA_forecast_data %>% 
  pull(.mean)

ARIMA_comparison <- tibble(
  Date = full_data %>% tail(8) %>% pull(Date),
  Actual = ARIMA_actual_values,
  Predicted = ARIMA_predicted_values
)

ARIMA_mse <- mean((ARIMA_comparison$Actual - ARIMA_comparison$Predicted)^2)
ARIMA_mae <- mean(abs(ARIMA_comparison$Actual - ARIMA_comparison$Predicted))
ARIMA_rmse <- sqrt(ARIMA_mse)

ARIMA_error_metrics <- tibble(
  MSE = ARIMA_mse,
  MAE = ARIMA_mae,
  RMSE = ARIMA_rmse
)
print(ARIMA_comparison)

## # A tibble: 8 × 3
##      Date Actual Predicted
##     <qtr>  <int>     <dbl>
## 1 2022 Q1   3509     3565.
## 2 2022 Q2   3432     3352.
## 3 2022 Q3   3357     3711.
## 4 2022 Q4   3850     4234.
## 5 2023 Q1   3328     3745.
## 6 2023 Q2   3290     3519.
## 7 2023 Q3   3225     3714.
## 8 2023 Q4   3677     4431.

print(ARIMA_error_metrics)

## # A tibble: 1 × 3
##       MSE   MAE  RMSE
##     <dbl> <dbl> <dbl>
## 1 164487.  345.  406.

best_ARIMA_model %>% 
  forecast(h = 8) %>%
  autoplot(full_data) + 
  labs(title = "ARIMA Retail Trade Forecasts", x = "Time", y = "Retail.Trade")

# Forecasting
forecast_NNAR <- fit_NNAR %>% forecast(h = 8)
# plot
autoplot(full_data) +
  autolayer(forecast_NNAR) +
  labs(title = "Retail Trade NNAR Forecasts", x = "Time", y = "Retail.Trade")

## Plot variable not specified, automatically selected `.vars = Retail.Trade`

accuracy(forecast_NNAR, full_data) %>% 
  select(.model, RMSE, MAE, MAPE, MASE)

## # A tibble: 1 × 5
##   .model                RMSE   MAE  MAPE  MASE
##   <chr>                <dbl> <dbl> <dbl> <dbl>
## 1 NNETAR(Retail.Trade)  389.  302.  8.84  3.29