The economic impact of international tourism in the United States is substantial, with global travelers contributing billions annually to the U.S. economy. Before the COVID-19 pandemic, international visitors spent approximately $233.5 billion in 2019 alone, which averages to nearly $640 million daily (International Trade Administration, n.d.). With U.S. tourism accounting for 14.5% of total global spending on international travel (Zaheer, 2023), understanding the trends and drivers behind these expenditures is vital for sustaining this economic sector.
While tourist spending is influenced by various socioeconomic factors, one of the most reliable indicators of such expenditures is the volume of international arrivals, commonly referred to as “Tourist Arrivals.” Fluctuations in arrivals often correlate closely with factors like currency exchange rates, seasonal trends, and global geopolitical stability (Hinton, 2024). In an attempt to explore the dynamics between tourist arrivals and spending, this project seeks to find the best predictive models of both metrics. Leveraging upon the function of time and SARIMA models, this study aims to provide a robust and comprehensive understanding of U.S. tourism’s economic contributions, as well as a tool for predicting future trends in international arrivals and spending.
Arrivals <- read.csv("Monthly_Tourist_Arrivals_R.csv")
Exports <- read.csv("Monthly_Exports_R.csv")
Arrivals$Total <- as.numeric(Arrivals$Total)
Arrivals$Time <- floor_date(as.Date(Arrivals$Time), "month")
Exports$Time <- as.Date(Exports$Time)
Masterdata0 <- merge(Arrivals, Exports, by = "Time", all = TRUE)For this project, we sourced two datasets from the International Trade Administration (ITA) website, both available in Excel format. The first dataset, “Monthly Exports Imports Balance,” was retrieved from the ITA’s International Travel Receipts and Payments Program. It provides monthly tourism-related trade data from January 1999 to July 2024, including metrics like total travel exports (receipts, or tourist spending) and imports (payments, or US citizen spending abroad), which are further subdivided into categories such as Travel Spending, Medical/ Education/ Workers Spending, and Passenger Fare Receipts. All figures are in millions of dollars and are seasonally adjusted.
The second dataset, “Monthly Arrivals 2000 to Present – Country of Residence (COR),” was acquired from the ITA’s I-94 Arrivals Program. This file contains monthly data on the number of tourists arriving in the U.S. from January 2000 to July 2024, with a breakdown by world region (e.g., Central America, Western Europe) and individual countries of residence. This dataset also includes year-over-year changes for each category. No specific filters or exclusions were applied to the data before downloading.
The datasets used in this study inherently exclude certain categories of visitors. For instance, the International Visitor Arrivals Program - ADIS I-94 does not count travelers who do not stay overnight in the U.S. or those who visit on certain visa types. This exclusion may introduce some problems, as short-term visitors, particularly those near borders or in airports, may still contribute significantly to tourism-related spending but are not included in official tourist arrival figures. Another potential bias is related to the seasonality of tourist arrivals, which can be influenced by holidays and business cycles. Significant anomalies like the COVID-19 pandemic have also impacted arrival and spending patterns.
To facilitate our analysis, we simplified both datasets by focusing on aggregated monthly data. For the “Monthly Exports Imports Balance” file, we retained only the total tourism-related exports (tourist spending). Similarly, for the “Monthly Arrivals” file, we retained only the total arrivals and subtotal information by world region, removing individual country data to streamline our analysis. The figure below shows a quick comparison of the two main time series remaining.
plot1 <- ggplot(Masterdata0, aes(x = Time, y = Total)) +
geom_line(color = "darkblue", size = 1) +
ggtitle("US Monthly Tourist Arrivals (2000-2024)") +
xlab("Month") +
ylab("Tourist Arrivals") +
scale_x_date(date_labels = "%b %Y", date_breaks = "36 months")+
theme_light()+
scale_y_continuous(labels = scales::label_number(scale_cut = cut_short_scale()))
plot2 <- ggplot(Masterdata0, aes(x = Time, y = Exports)) +
geom_line(color = "darkred", size = 1) +
ggtitle("US Monthly Tourist Exports, seasonally adjusted (2000-2024)") +
xlab("Month") +
ylab("Export Values (million $)") +
scale_x_date(date_labels = "%b %Y", date_breaks = "36 months")+
theme_light()+
scale_y_continuous(labels = scales::label_number(scale_cut = cut_short_scale()))
grid.arrange(plot1, plot2, nrow=2)
Given the unique impact of COVID-19 on travel and spending, we intentionally excluded data post-January 2020. This choice allows us to analyze trends under more typical conditions, thus providing a more stable basis for examining seasonal and trend-based patterns in tourist arrivals and spending from 2000 to 2019. The figure below illustrates such exclusion of data.
Masterdata <-Masterdata0 %>%
filter(Time < as.Date("2020-01-01"))
Arrivals_ts = ts(Masterdata$Total,start= c(2000,1), frequency = 12)
Exports_ts = ts(Masterdata$Exports,start= c(2000,1), frequency = 12)
tstools::tsplot(list("Arrivals (millions), left axis"=Arrivals_ts/1000000),
tsr = list("Exports (billions $), right axis" =Exports_ts/1000),
theme = init_tsplot_theme(line_colors = c("darkblue", "darkred")),
plot_title = "US Monthly Tourist Arrivals & Spending (2000-2019)")
To conduct a rigorous statistical analysis of U.S. monthly tourist arrivals and spending, we assume that each account in our data is independent of one another and meets all criteria described by the International Trade Administration for their corresponding reporting data. The following paragraphs outline our general approach to finding the best models for our data.
Before fitting any model, we will visualize the data with a time series plot and compute the Autocorrelation Function (ACF) and Partial Autocorrelation Function (PACF) plots to understand the underlying patterns if appropriate. We will then perform several stationarity tests, including the Augmented Dickey-Fuller (ADF) test, Phillips-Perron (PP) test, Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test, and the Elliott-Rothenberg-Stock (ERS) test. More information on these tests and beyond will be available in the Appendix.
After fitting the models, we will validate their effectiveness by plotting the fitted values with the original data and then analyzing the residuals. We will check for independence, normality, and stationarity of the residuals, employing tests such as the Shapiro-Wilk, Anderson-Darling, Jarque-Bera, and Lilliefors tests to ensure they follow a normal distribution. Additionally, ACF and density plots, Box-Ljung test, and other stationarity tests will inform us about potential issues with correlation and stationarity. In general, if the residuals are stationary, independent, and normally distributed, the model’s predictions can be deemed reliable.
When comparing models, we will utilize the Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC), and R-squared values (for linear models). A lower AIC and BIC, along with a higher R-squared value, indicate a better-fitting model. If residuals are stationary, we may enhance the model by incorporating additional terms to improve its predictive capability.
For validation, we will implement a cross-validation approach by excluding data from 2019, training our models with the available data up to 2018, and subsequently comparing the model’s predictions against actual arrivals and spending for that year. Finally, we will forecast tourist arrivals and spending for two years into the future, assuming the absence of COVID-19-related disruptions.
In the event of identifying non-stationarity due to trends within the time series, we will apply function of time models to detrend the data. Specifically, we will utilize regression techniques to fit linear, quadratic, sinusoidal, and periodic trends with seasonal means. By identifying the appropriate model that best captures these trends, we will aim to derive residuals that should exhibit stationary behavior—constant variance and no serial correlation. Within this regard, we will also utilize decomposition techniques and plot accordingly to compare with our function of time models.
Seasonal Autoregressive Integrated Moving Average (SARIMA) models will be employed to account for seasonal behavior in the data. SARIMA models are denoted as SARIMA(p, d, q) × (P, D, Q) [s], where: p: Order of the autoregressive component d: Degree of differencing q: Order of the moving average component P: Seasonal autoregressive order D: Seasonal differencing order Q: Seasonal moving average order s: Length of the seasonal cycle
Prior to fitting a SARIMA model, we will first difference the data (normally and/or seasonally) as necessary to ensure stationarity. Different models will be fitted (based on ACF and PACF plots), validated by series of statistical tests, and analyzed in terms of residuals to confirm that they are stationary and uncorrelated. If the residuals follow a normal distribution, we can confidently utilize the model for forecasting future values.
To facilitate this analysis, I have self-defined new functions in R to streamline various processes, including:
check_stationarity(): This function performs stationarity tests (ADF, PP, KPSS, ERS) and summarizes the results in a user-friendly table.
# Stationarity Tests
check_stationarity <- function(data) {
adf_result <- adf.test(data)
pp_result <- pp.test(data)
kpss_result <- kpss.test(data)
urers_result <- ur.ers(data, type = "DF-GLS", model = "constant")
results <- data.frame(
Test = c("ADF Test", "Phillips-Perron Test", "KPSS Test", "ERS Test"),
`Test Statistic` = c(round(adf_result[["statistic"]][["Dickey-Fuller"]], 3),
round(pp_result[["statistic"]][["Dickey-Fuller Z(alpha)"]], 3),
round(kpss_result[["statistic"]][["KPSS Level"]], 3),
round(urers_result@teststat, 3)),
`P-Value` = c(round(adf_result[["p.value"]],2),
round(pp_result[["p.value"]],2),
round(kpss_result[["p.value"]],2),
"N/A"),
Interpretation = c(
ifelse(adf_result[["p.value"]] < 0.05, "Stationary", "Non-Stationary"),
ifelse(pp_result[["p.value"]] < 0.05, "Stationary", "Non-Stationary"),
ifelse(kpss_result[["p.value"]] < 0.05, "Non-Stationary", "Stationary"),
ifelse(urers_result@teststat < urers_result@cval[2], "Stationary", "Non-Stationary")))
kable(results, align = 'c') %>%
kable_styling(full_width = F, bootstrap_options = c("striped", "hover", "condensed")) %>%
add_header_above(c("Stationarity Test Summary" = 4))
}four_plots() and two_plots(): These functions generate ACF and PACF plots for different segments of the data to visually assess autocorrelation.
# ACF/PACF/First-Last Third Plotting
four_plots <- function(data) {
n <- length(data)
first_third <- data[1:floor(n / 3)]
last_third <- data[floor(2 * n / 3):n]
par(mfrow = c(2, 2))
acf(data, main = "ACF of Whole Data", lag.max =30)
pacf(data, main = "PACF of Whole Data", lag.max =30)
acf(first_third, main = "ACF of First Third")
acf(last_third, main = "ACF of Last Third")
par(mfrow = c(1, 1))
}
# ACF/PACF Plotting
two_plots <- function(data) {
par(mfrow = c(1, 2))
acf(data, main = "ACF of Whole Data", lag.max =50)
pacf(data, main = "PACF of Whole Data", lag.max =50)
par(mfrow = c(1, 1))
}
# Differencing Plotting
diff_plots <- function(data, name) {
layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
astsa::tsplot(data, main= name)
abline(h = mean(data, na.rm = TRUE), col = "red", lty = 2)
acf(data, lag.max =30)
pacf(data, lag.max =30)
par(mfrow = c(1,1))
}model_summary_table() and sarima_summary_table( ): These functions summarize the key metrics of fitted models, including R^2 , AIC and BIC values, to facilitate comparison.
# Model Summary
model_summary_table <- function(models_list) {
model_names <- names(models_list)
if (is.null(model_names)) {
model_names <- paste("Model", seq_along(models_list)) }
summary_df <- data.frame(Model = character(),
`Adjusted $R^2$ ` = numeric(),
AIC = numeric(),
BIC = numeric(),
stringsAsFactors = FALSE)
for (i in seq_along(models_list)) {
model <- models_list[[i]]
model_summary <- summary(model)
summary_df[i, ] <- c(
model_names[i],
round(model_summary[["adj.r.squared"]], 2),
round(AIC(model), 2),
round(BIC(model), 2)) }
kable(summary_df, align = 'c', col.names = c("Model", "Adjusted R²", "AIC", "BIC")) %>%
kable_styling(full_width = F, bootstrap_options = c("striped", "hover", "condensed")) %>%
add_header_above(c("Model Summary" = 4))
}
# SARIMA Model Summary
sarima_summary_table <- function(models_list) {
model_names <- names(models_list)
if (is.null(model_names)) {
model_names <- paste("Model", seq_along(models_list)) }
summary_df <- data.frame(Model = character(),
AIC = numeric(),
BIC = numeric(),
stringsAsFactors = FALSE)
for (i in seq_along(models_list)) {
model <- models_list[[i]]
model_summary <- summary(model)
summary_df[i, ] <- c(
model_names[i],
round(AIC(model), 2),
round(BIC(model), 2)) }
kable(summary_df, align = 'c', col.names = c("Model", "AIC", "BIC")) %>%
kable_styling(full_width = F, bootstrap_options = c("striped", "hover", "condensed")) %>%
add_header_above(c("Model Summary" = 3))
}check_normality_residuals(): This function evaluates the residuals from fitted models through a series of different tests to ensure they conform to normality assumptions.
# Normality Tests Function
check_normality_residuals <- function(model) {
residuals <- residuals(model)
shapiro_result <- shapiro.test(residuals)
ad_result <- ad.test(residuals)
jb_result <- jarque.bera.test(residuals)
lillie_result <- lillie.test(residuals)
skew <- skewness(residuals)
kurt <- kurtosis(residuals)
# Summary table
results <- data.frame(
Test = c("Shapiro-Wilk Test", "Anderson-Darling Test", "Jarque-Bera Test", "Lilliefors Test", "Skewness", "Kurtosis"),
`Test Statistic` = c(round(shapiro_result$statistic, 3),
round(ad_result$statistic, 3),
round(jb_result$statistic, 3),
round(lillie_result$statistic, 3),
round(skew, 2),
round(kurt, 2)),
`P-Value` = c(round(shapiro_result$p.value, 3),
round(ad_result$p.value, 3),
round(jb_result$p.value, 3),
round(lillie_result$p.value, 3),
"N/A", "N/A"),
Interpretation = c(
ifelse(shapiro_result$p.value > 0.05, "Normal", "Non-Normal"),
ifelse(ad_result$p.value > 0.05, "Normal", "Non-Normal"),
ifelse(jb_result$p.value > 0.05, "Normal", "Non-Normal"),
ifelse(lillie_result$p.value > 0.05, "Normal", "Non-Normal"),
ifelse(abs(skew) < 0.5, "Approximately Symmetric", ifelse(skew > 0, "Right-Skewed", "Left-Skewed")),
ifelse(abs(kurt - 3) < 0.5, "Mesokurtic", ifelse(kurt > 3, "Leptokurtic (Heavy Tails)", "Platykurtic (Light Tails)"))
)
)
kable(results, align = 'c') %>%
kable_styling(full_width = F, bootstrap_options = c("striped", "hover", "condensed")) %>%
add_header_above(c("Normality Test Summary" = 4))
} # Time Series Summary Table
ts_summary_table <- function(ts_list) {
ts_names <- names(ts_list)
if (is.null(ts_names)) {
ts_names <- paste("Series", seq_along(ts_list)) }
summary_df <- data.frame(Series = character(),
Start = character(),
End = character(),
Min = numeric(),
Max = numeric(),
Mean = numeric(),
Median = numeric(),
stringsAsFactors = FALSE)
for (i in seq_along(ts_list)) {
ts_data <- ts_list[[i]]
summary_df[i, ] <- c(
ts_names[i],
start(ts_data) %>% paste(collapse = "-"),
end(ts_data) %>% paste(collapse = "-"),
round(min(ts_data, na.rm = TRUE), 2),
round(max(ts_data, na.rm = TRUE), 2),
round(mean(ts_data, na.rm = TRUE), 2),
round(median(ts_data, na.rm = TRUE), 2)
)
}
kable(summary_df, align = 'c',
col.names = c("Series", "Start", "End", "Min", "Max", "Mean", "Median")) %>%
kable_styling(full_width = F, bootstrap_options = c("striped", "hover", "condensed")) %>%
add_header_above(c("Time Series Summary" = 7))
}calculate_r2 <- function(predicted_values, actual_values) {
ss_residual <- sum((actual_values - predicted_values)^2)
ss_total <- sum((actual_values - mean(actual_values))^2)
r_squared <- 1 - (ss_residual / ss_total)
cat("R-squared of the model is:", round(r_squared, 4), "\n")}These functions enhance the analysis process by providing clarity and efficiency. Each function’s details and the tests used are elaborated upon in the appendix.
U.S. monthly tourist arrivals showed a notable increase over the observed period, beginning with 2,866,229 visitors in January 2000 and reaching 6,899,661 by July 2024, indicating an overall upward trend. The plot below illustrates the monthly arrival data from 2000 to July 2024, segmented by different world regions as well as the total aggregate. The mean monthly arrivals were approximately 4,709,139, with a median of 4,817,910. Notably, Western Europe accounted for an average of 902,992 monthly arrivals, Asia for 618,172, and South America for 301,831.
The lowest recorded arrivals occurred in April 2020, when the number of visitors dropped sharply to 248,486, corresponding to the onset of U.S. border closures and social distancing measures due to the COVID-19 pandemic (Gonzalez-Barrera, 2020). Conversely, the highest number of monthly arrivals was recorded in August 2014, with 8,418,370 visitors.
astsa::tsplot(x= Masterdata0$Time, y= Masterdata0$Total/1000000,
col = "darkblue",
lwd = 2,
main = "US Monthly Tourist Arrivals (2000-2024)",
xlab = "Year",
ylab = "Tourist Arrivals (millions)")
lines(x= Masterdata0$Time, y= Masterdata0$Africa/1000000, col = 1, lwd = 1.5)
lines(x= Masterdata0$Time, y= Masterdata0$Asia/1000000, col = 2, lwd = 1.5)
lines(x= Masterdata0$Time, y= Masterdata0$W.Europe/1000000, col = 3, lwd = 1.5)
lines(x= Masterdata0$Time, y= Masterdata0$E.Europe/1000000, col = 4, lwd = 1.5)
lines(x= Masterdata0$Time, y= Masterdata0$C.America/1000000, col = 5, lwd = 1.5)
lines(x= Masterdata0$Time, y= Masterdata0$Caribbean/1000000, col = 6, lwd = 1.5)
lines(x= Masterdata0$Time, y= Masterdata0$S.America/1000000, col = 7, lwd = 1.5)
lines(x= Masterdata0$Time, y= Masterdata0$Oceania/1000000, col = 8, lwd = 1.5)
lines(x= Masterdata0$Time, y= Masterdata0$M.East/1000000, col = 9, lwd = 1.5)
legend("topleft", col=c("darkblue",1:9), lty=1, lwd=2, pch=20,
legend=c("Total", "Africa", "Asia", "W.Europe", "E.Europe", "C.America", "Caribbean", "S.America", "Oceania","M.East"), bg="white",ncol=2)
Seasonal patterns are evident in the plot above, with clear peaks and troughs repeating annually. The plot below signifies that the summer months of July and August consistently show the highest arrivals, while January and February experience the lowest. This trend aligns with travel behavior, as summer weather is more favorable for tourism, and winter months see a decrease in travel as people return to work or school after the holidays. Additionally, colder weather and snow in northern regions may deter travel during these months.
Masterdata <- Masterdata %>%
select("Time","Total","Exports") %>%
mutate("Month" = month(Time)) %>%
mutate("Year" = year(Time))
p <- ggplot(Masterdata, aes(x=Month, y=Total/1000000, color = Year)) +
geom_line(alpha = 0.7) +
labs(x = "Month",
y = "Total Arrivals (millions)",
title = "Total Arrivals by Month") +
theme_minimal() +
scale_x_continuous(breaks=c(2,4, 6, 8,10,12)) +
transition_states(Year)
animate(p, renderer = gifski_renderer())
Monthly tourist spending, or “Tourist Exports,” also shows a marked increase, starting at 9.2 billion dollars in January 2000 and reaching 21.1 billion dollars in July 2024. The mean monthly spending over this period was 13.7 billion dollars, with a median of 13.1 billion dollars. The lowest recorded spending occurred in September 2020, at 3.7 billion dollars, while the highest was in May 2024, with 21.66 billion dollars. However, it is noteworthy that this dataset is seasonally adjusted, which means that nominal spending figures may vary; for instance, the lowest nominal spending month may have occurred earlier in 2020 when the tourism sector was effectively shut down.
For additional context on seasonal adjustments, please refer to the Appendix, where the seasonal adjustment methodology is outlined.
ggplot(Masterdata, aes(x = Time, y = Exports)) +
geom_line(color = "darkred", size = 1) +
ggtitle("US Monthly Tourist Exports, seasonally adjusted (2000-2024)") +
xlab("Month") +
ylab("Export Values (million $)") +
scale_x_date(date_labels = "%b %Y", date_breaks = "36 months")+
theme_light()+
scale_y_continuous(labels = scales::label_number(scale_cut = cut_short_scale()))
The table below compares statistical summaries of monthly tourist arrivals and exports for the entire period (2000-2024) versus the period excluding data from January 2020 onward.
ar_full <- ts(Masterdata0$Total,start= c(2000,1), frequency = 12)
ex_full <- ts(Masterdata0$Exports,start= c(2000,1), frequency = 12)
ts_series <- list("Tourist Arrivals (2000-2024) (millions)" = ar_full/1000000,
"Tourist Arrivals (2000-2019) (millions)" = Arrivals_ts/1000000,
"Tourist Spending (2000-2024) (billions USD)" = ex_full/1000,
"Tourist Spending (2000-2019) (billions USD)" = Exports_ts/1000)
ts_summary_table(ts_series)| Series | Start | End | Min | Max | Mean | Median |
|---|---|---|---|---|---|---|
| Tourist Arrivals (2000-2024) (millions) | 2000-1 | 2024-7 | 0.25 | 8.42 | 4.71 | 4.82 |
| Tourist Arrivals (2000-2019) (millions) | 2000-1 | 2019-12 | 2.09 | 8.42 | 4.96 | 4.91 |
| Tourist Spending (2000-2024) (billions USD) | 2000-1 | 2024-7 | 3.72 | 21.66 | 13.7 | 13.11 |
| Tourist Spending (2000-2019) (billions USD) | 2000-1 | 2019-12 | 6.7 | 20.82 | 13.86 | 12.89 |
In total, we evaluated five different models to capture trends in the total arrivals time series, incorporating linear, quadratic, seasonal, and sinusoidal components. In order to mitigate the multicollinearity problems associated with VIF, we intentionally used the centered time object in all of our models, if appropriate.
In this model, we centered the time variable and fit a linear trend to arrivals. The output reveals a significant upward trend, with an estimated increase of 215,025 arrivals per time unit (p < 0.05) and an intercept of 4,958,137. Residuals are centered and show a balanced spread. An adjusted R-squared of 0.7164 indicates that 71.64% of the variance in arrivals is explained by time, and a highly significant F-statistic (604.7, p < 0.05) supports the model’s strength and the importance of time in explaining arrival trends. While the residuals appear stationary, they are not independent and do not follow a normal distribution.
t <- time(Arrivals_ts)
t_centered <- t - mean(t)
# Linear
ar_lm <- lm(Arrivals_ts ~ t_centered)
summary(ar_lm)
Call:
lm(formula = Arrivals_ts ~ t_centered)
Residuals:
Min 1Q Median 3Q Max
-1828757 -525528 -55961 354770 2465742
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4958137 50483 98.22 <2e-16 ***
t_centered 215025 8744 24.59 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 782100 on 238 degrees of freedom
Multiple R-squared: 0.7176, Adjusted R-squared: 0.7164
F-statistic: 604.7 on 1 and 238 DF, p-value: < 2.2e-16ts.plot(cbind(Arrivals_ts/1000000,ar_lm$fitted.values/1000000), col=c("darkblue", "darkred"),
ylab = "Arrivals (millions)", main = "Original Data with Fitted Values - Linear", lwd =1.5)
legend("topleft", legend=c("Original Data", "Fitted Values"), col=c("darkblue", "darkred"), lty=1:1,lwd=2)
checkresiduals(ar_lm$residuals)
Ljung-Box test
data: Residuals
Q* = 176.1, df = 10, p-value < 2.2e-16
Model df: 0. Total lags used: 10check_stationarity(ar_lm$residuals)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| ADF Test | -5.456 | 0.01 | Stationary |
| Phillips-Perron Test | -113.997 | 0.01 | Stationary |
| KPSS Test | 0.109 | 0.1 | Stationary |
| ERS Test | -7.740 | N/A | Stationary |
check_normality_residuals(ar_lm)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| Shapiro-Wilk Test | 0.973 | 0 | Non-Normal |
| Anderson-Darling Test | 2.066 | 0 | Non-Normal |
| Jarque-Bera Test | 13.541 | 0.001 | Non-Normal |
| Lilliefors Test | 0.084 | 0 | Non-Normal |
| Skewness | 0.570 | N/A | Right-Skewed |
| Kurtosis | 3.250 | N/A | Mesokurtic |
In this model, we included both a linear term (t_centered) and a quadratic term (I(t^2)) to fit arrivals. The output indicates that neither the intercept, linear, nor quadratic coefficients are statistically significant (p > 0.05), suggesting that the quadratic term does not add explanatory power to the model. However, an adjusted R-squared at 0.7152 means 71.52% of the variance in arrivals is still explained by the model. This is lower than the linear model alone, so we will not keep the quadratic term in following models. Meanwhile, a highly significant F-statistic (301.1, p < 0.05) confirms that the model overall is meaningful, despite the individual terms lacking significance. Residuals are stationary but lack independence and normality.
# Quadratic
ar_quad <- lm(Arrivals_ts ~ t_centered + I(t^2) )
summary(ar_quad)
Call:
lm(formula = Arrivals_ts ~ t_centered + I(t^2))
Residuals:
Min 1Q Median 3Q Max
-1818018 -527921 -57131 355228 2463173
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.738e+08 6.855e+09 0.127 0.899
t_centered 1.080e+06 6.821e+06 0.158 0.874
I(t^2) -2.151e+02 1.697e+03 -0.127 0.899
Residual standard error: 783700 on 237 degrees of freedom
Multiple R-squared: 0.7176, Adjusted R-squared: 0.7152
F-statistic: 301.1 on 2 and 237 DF, p-value: < 2.2e-16ts.plot(cbind(Arrivals_ts/1000000,ar_quad$fitted.values/1000000), col=c("darkblue", "darkred"),
ylab = "Arrivals (millions)", main = "Original Data with Fitted Values - Quadratic", lwd =1.5)
legend("topleft", legend=c("Original Data", "Fitted Values"), col=c("darkblue", "darkred"), lty=1:1,lwd=2)
checkresiduals(ar_quad$residuals)
Ljung-Box test
data: Residuals
Q* = 176.18, df = 10, p-value < 2.2e-16
Model df: 0. Total lags used: 10check_stationarity(ar_quad$residuals)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| ADF Test | -5.480 | 0.01 | Stationary |
| Phillips-Perron Test | -113.988 | 0.01 | Stationary |
| KPSS Test | 0.107 | 0.1 | Stationary |
| ERS Test | -7.724 | N/A | Stationary |
check_normality_residuals(ar_quad)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| Shapiro-Wilk Test | 0.973 | 0 | Non-Normal |
| Anderson-Darling Test | 2.087 | 0 | Non-Normal |
| Jarque-Bera Test | 13.648 | 0.001 | Non-Normal |
| Lilliefors Test | 0.083 | 0 | Non-Normal |
| Skewness | 0.570 | N/A | Right-Skewed |
| Kurtosis | 3.240 | N/A | Mesokurtic |
In this model, we included a centered time trend (t_centered) and monthly indicator variables (t_month) to capture both linear trends and seasonal patterns in arrivals. The output shows that both the intercept and time trend are highly significant (p < 0.05), indicating a clear upward trend over time. Many monthly coefficients are also significant, suggesting distinct seasonal effects on arrivals (e.g., increased arrivals in July and August). The adjusted R-squared is 0.9131, meaning 91.31% of the variance in arrivals is explained by the model, a marked improvement in explanatory power. A highly significant F-statistic (210.4, p < 0.05) confirms the model’s strong overall fit. However, the residuals show non-stationarity, non-independence, and non-normality.
# Linear + Seasonal Trend (0.9128 with t^2 < 0.9131 without t^2)
t_month <- factor(month(Masterdata$Time))
ar_month <- lm(Arrivals_ts ~ t_centered + t_month)
summary(ar_month)
Call:
lm(formula = Arrivals_ts ~ t_centered + t_month)
Residuals:
Min 1Q Median 3Q Max
-1123838 -251167 14405 276126 1177669
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4279501 96805 44.207 < 2e-16 ***
t_centered 213315 4845 44.027 < 2e-16 ***
t_month2 -402561 136868 -2.941 0.00361 **
t_month3 594818 136870 4.346 2.09e-05 ***
t_month4 845284 136873 6.176 3.01e-09 ***
t_month5 710543 136877 5.191 4.64e-07 ***
t_month6 451356 136882 3.297 0.00113 **
t_month7 1854503 136889 13.547 < 2e-16 ***
t_month8 2018620 136897 14.746 < 2e-16 ***
t_month9 678724 136906 4.958 1.40e-06 ***
t_month10 622646 136916 4.548 8.83e-06 ***
t_month11 181691 136927 1.327 0.18587
t_month12 587998 136940 4.294 2.60e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 432800 on 227 degrees of freedom
Multiple R-squared: 0.9175, Adjusted R-squared: 0.9131
F-statistic: 210.4 on 12 and 227 DF, p-value: < 2.2e-16ts.plot(cbind(Arrivals_ts/1000000,ar_month$fitted.values/1000000 ), col=c("darkblue", "darkred"),
ylab = "Arrivals (millions)", main = "Original Data with Fitted Values - Linear + Seasonal", lwd =1.5)
legend("topleft", legend=c("Original Data", "Fitted Values"), col=c("darkblue", "darkred"), lty=1:1,lwd=2)
checkresiduals(ar_month$residuals)
Ljung-Box test
data: Residuals
Q* = 713.33, df = 10, p-value < 2.2e-16
Model df: 0. Total lags used: 10check_stationarity(ar_month$residuals)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| ADF Test | -3.270 | 0.08 | Non-Stationary |
| Phillips-Perron Test | -44.256 | 0.01 | Stationary |
| KPSS Test | 0.233 | 0.1 | Stationary |
| ERS Test | -1.551 | N/A | Non-Stationary |
check_normality_residuals(ar_month)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| Shapiro-Wilk Test | 0.992 | 0.204 | Normal |
| Anderson-Darling Test | 0.757 | 0.048 | Non-Normal |
| Jarque-Bera Test | 0.225 | 0.894 | Normal |
| Lilliefors Test | 0.063 | 0.022 | Non-Normal |
| Skewness | 0.010 | N/A | Approximately Symmetric |
| Kurtosis | 3.150 | N/A | Mesokurtic |
In this model, we fit a linear trend (t_centered) alongside a seasonal cosine term (cosine_t) to capture both the trend and periodic seasonal effects in arrivals. The intercept and time trend are highly significant (p < 0.05), suggesting a consistent upward trend, with an estimated increase of 215,116 arrivals per time unit. However, the cosine term is not significant (p = 0.155), indicating it does not add meaningful explanatory power for seasonality. The adjusted R-squared value of 0.7176 shows that the model explains 71.76% of the variance, while a highly significant F-statistic (304.7, p < 0.05) confirms the model’s overall fit. Residuals are stationary, yet they lack independence and normality.
# Linear + Cosine Trend
cosine_t <- cos(2 * pi * t / 12)
ar_cos <- lm(Arrivals_ts ~ t_centered + cosine_t)
summary(ar_cos)
Call:
lm(formula = Arrivals_ts ~ t_centered + cosine_t)
Residuals:
Min 1Q Median 3Q Max
-1852455 -526048 -64173 376932 2527897
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4976231 51946 95.797 <2e-16 ***
t_centered 215116 8725 24.655 <2e-16 ***
cosine_t -109418 76702 -1.427 0.155
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 780400 on 237 degrees of freedom
Multiple R-squared: 0.72, Adjusted R-squared: 0.7176
F-statistic: 304.7 on 2 and 237 DF, p-value: < 2.2e-16ts.plot(cbind(Arrivals_ts/1000000,ar_cos$fitted.values/1000000 ), col=c("darkblue", "darkred"),
ylab = "Arrivals (millions)", main = "Original Data with Fitted Values - Linear + Cosine", lwd =1.5)
legend("topleft", legend=c("Original Data", "Fitted Values"), col=c("darkblue", "darkred"), lty=1:1,lwd=2)
checkresiduals(ar_cos$residuals)
Ljung-Box test
data: Residuals
Q* = 177.01, df = 10, p-value < 2.2e-16
Model df: 0. Total lags used: 10check_stationarity(ar_cos$residuals)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| ADF Test | -5.472 | 0.01 | Stationary |
| Phillips-Perron Test | -114.949 | 0.01 | Stationary |
| KPSS Test | 0.103 | 0.1 | Stationary |
| ERS Test | -7.842 | N/A | Stationary |
check_normality_residuals(ar_cos)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| Shapiro-Wilk Test | 0.975 | 0 | Non-Normal |
| Anderson-Darling Test | 1.907 | 0 | Non-Normal |
| Jarque-Bera Test | 14.711 | 0.001 | Non-Normal |
| Lilliefors Test | 0.078 | 0.001 | Non-Normal |
| Skewness | 0.590 | N/A | Right-Skewed |
| Kurtosis | 3.310 | N/A | Mesokurtic |
In this model, we include both a linear time trend (t_centered) and seasonal components (cosine_t and t_month) to capture the trend and monthly variations in tourist arrivals. The intercept and trend term are highly significant (p < 0.05), with an estimated increase of 213,407 arrivals per time unit. The cosine term and many monthly indicators are also significant, suggesting both trend and seasonal effects are present in the data. The high adjusted R-squared of 0.9154 indicates the model explains 91.54% of the variance in arrivals, supported by the significant F-statistic (199.9, p < 0.05) for overall model fit. Residuals are largely normal and stationary, though they show some autocorrelation, indicating residual dependence over time.
# Linear + Seasonal + Cosine
ar_both <- lm(Arrivals_ts ~ t_centered + cosine_t + t_month)
summary(ar_both)
Call:
lm(formula = Arrivals_ts ~ t_centered + cosine_t + t_month)
Residuals:
Min 1Q Median 3Q Max
-1035888 -297866 9284 274209 1196632
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4297509 95791 44.864 < 2e-16 ***
t_centered 213407 4782 44.625 < 2e-16 ***
cosine_t -111170 41988 -2.648 0.008676 **
t_month2 -402376 135092 -2.979 0.003212 **
t_month3 595153 135093 4.405 1.63e-05 ***
t_month4 845736 135096 6.260 1.91e-09 ***
t_month5 711075 135100 5.263 3.28e-07 ***
t_month6 451933 135106 3.345 0.000963 ***
t_month7 1855090 135112 13.730 < 2e-16 ***
t_month8 2019183 135120 14.944 < 2e-16 ***
t_month9 679226 135129 5.027 1.02e-06 ***
t_month10 623052 135139 4.610 6.72e-06 ***
t_month11 181965 135150 1.346 0.179523
t_month12 588107 135162 4.351 2.05e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 427200 on 226 degrees of freedom
Multiple R-squared: 0.92, Adjusted R-squared: 0.9154
F-statistic: 199.9 on 13 and 226 DF, p-value: < 2.2e-16ts.plot(cbind(Arrivals_ts/1000000,ar_both$fitted.values/1000000 ), col=c("darkblue", "darkred"),
ylab = "Arrivals (millions)", main = "Original Data with Fitted Values - Linear + Seasonal + Cosine", lwd =1.5)
legend("topleft", legend=c("Original Data", "Fitted Values"), col=c("darkblue", "darkred"), lty=1:1,lwd=2)
checkresiduals(ar_both$residuals)
Ljung-Box test
data: Residuals
Q* = 694.36, df = 10, p-value < 2.2e-16
Model df: 0. Total lags used: 10check_stationarity(ar_both$residuals)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| ADF Test | -3.193 | 0.09 | Non-Stationary |
| Phillips-Perron Test | -45.899 | 0.01 | Stationary |
| KPSS Test | 0.226 | 0.1 | Stationary |
| ERS Test | -1.676 | N/A | Non-Stationary |
check_normality_residuals(ar_both)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| Shapiro-Wilk Test | 0.996 | 0.805 | Normal |
| Anderson-Darling Test | 0.229 | 0.809 | Normal |
| Jarque-Bera Test | 0.765 | 0.682 | Normal |
| Lilliefors Test | 0.035 | 0.665 | Normal |
| Skewness | 0.130 | N/A | Approximately Symmetric |
| Kurtosis | 2.890 | N/A | Mesokurtic |
Of these, two models with seasonal adjustments, “Linear + Seasonal” and “Linear + Seasonal + Cosine,” demonstrated the best performance in terms of trend accuracy and seasonal variability. These models also achieved the lowest AIC values (6925.19 and 6919.86), lowest BIC values (6973.92 and 6972.07), and the highest R-squared values (0.91 and 0.92). Based on these criteria, the “Linear + Seasonal + Cosine” model is the preferred fit for the tourist arrivals data.
# Summary
ar_models <- list("Linear Model" = ar_lm,
"Quadratic Model" = ar_quad,
"Linear + Seasonal" = ar_month,
"Linear + Cosine" = ar_cos,
"Linear + Seasonal + Cosine" = ar_both)
model_summary_table(ar_models)| Model | Adjusted R² | AIC | BIC |
|---|---|---|---|
| Linear Model | 0.72 | 7198.54 | 7208.98 |
| Quadratic Model | 0.72 | 7200.52 | 7214.45 |
| Linear + Seasonal | 0.91 | 6925.19 | 6973.92 |
| Linear + Cosine | 0.72 | 7198.49 | 7212.41 |
| Linear + Seasonal + Cosine | 0.92 | 6919.86 | 6972.07 |
However, it is important to note that this model’s residuals do not meet assumptions of stationarity, normality, or independence, which limits our confidence in the forecast intervals.
# Check model
check_model(ar_month)
Here, we try to decompose the time series. As shown in the plots, there are significant seasonality of frequency 12 months.
# Classical Seasonal Decomposition by Moving Averages
plot(decompose(Arrivals_ts/1000000))
# STL Decomposition
plot(stl(Arrivals_ts/1000000, s.window = 12))
Here, we trained the model with data up to 2018, and then tested against actual data in 2019. This technique is called Leave One Out Cross Validation (LOOCV) (Olivo, 2024). The predicted values and intervals miss the actual data almost every time and the \(R^2\) is negative, suggesting a poor fit.
# Train on Data Up to 2018 and Predict for 2019 (with 95% CI)
train_data <- window(Arrivals_ts, end = c(2018, 12))
test_data <- window(Arrivals_ts, start = c(2019, 1), end = c(2019, 12))
t_train <- time(train_data)
t_month_train <- t_month[1:228]
cosine_t_train <- cosine_t[1:228]
model_train <- lm(train_data ~ t_train + cosine_t_train + t_month_train)
t_test <- time(test_data)
t_month_test <- t_month[229:240]
cosine_t_test <- cosine_t[229:240]
predictions_2019 <- predict(model_train,
newdata = data.frame(t_train = t_test,
cosine_t_train = cosine_t_test,
t_month_train = t_month_test),
interval = "confidence", level = 0.95)
# Combine actual and predicted data
time_full_2019 <- data.frame(
Date = c(time(train_data), time(test_data)),
Actual = c(as.numeric(train_data), as.numeric(test_data)),
Predicted = c(rep(NA, length(train_data)), predictions_2019[,"fit"]),
Lower_CI = c(rep(NA, length(train_data)), predictions_2019[,"lwr"]),
Upper_CI = c(rep(NA, length(train_data)), predictions_2019[,"upr"]))
time_only_2019 <- data.frame(
Date = as.Date(time(test_data)),
Actual = test_data,
Predicted = predictions_2019[,"fit"],
Lower_CI = predictions_2019[,"lwr"],
Upper_CI = predictions_2019[,"upr"])
ar_pred_2019 <- ggplot(time_only_2019, aes(x = Date)) +
geom_line(aes(y = Actual/1000000, color = "Actual")) +
geom_line(aes(y = Predicted/1000000, color = "Predicted")) +
geom_ribbon(aes(ymin = Lower_CI/1000000, ymax = Upper_CI/1000000), fill = "grey30", alpha = 0.4) +
labs(title = "Actual and Predicted Arrivals (2019) with 95% CI",
y = "Tourist Arrivals (millions)")+
scale_x_date(date_labels = "%b %Y", date_breaks = "2 months") +
theme_minimal() +
scale_color_manual(values = c("Actual" = "darkblue", "Predicted" = "red"))
ar_pred_full <- ggplot(time_full_2019, aes(x = Date)) +
geom_line(aes(y = Actual/1000000, color = "Actual")) +
geom_line(aes(y = Predicted/1000000, color = "Predicted")) +
geom_ribbon(aes(ymin = Lower_CI/1000000, ymax = Upper_CI/1000000), fill = "grey30", alpha = 0.4) +
labs(title = "Actual and Predicted Arrivals (2019) with 95% CI",
y = "Tourist Arrivals (millions)") +
theme_minimal() +
scale_color_manual(values = c("Actual" = "darkblue", "Predicted" = "red"))
grid.arrange(ar_pred_full, ar_pred_2019, nrow=2)
# R^2
calculate_r2(as.numeric(predictions_2019),as.numeric(test_data))R-squared of the model is: -0.675 Based on the poor result in the cross-validation test, when we applied the model to forecast arrivals for hypothetical years 2020 and 2021 (assuming no impact from COVID-19), the predictions likely exceeded what actual values would have been in such a scenario.
# Fit Model on Full Data and Forecast Two Years Ahead (with 95% CI)
ar_pred <- lm(Arrivals_ts ~ t + cosine_t + t_month)
future_dates <- ts(rep(NA, 24),start= c(2020,1), frequency = 12)
t_future <- time(future_dates)
cosine_t_future <- cos(2 * pi * time(future_dates) / 12)
t_month_future <- factor(month(seq(as.Date("2020-01-01"),
as.Date("2021-12-01"), by = "month")))
predictions_future <- predict(ar_pred,
newdata = data.frame(t = t_future,
cosine_t = cosine_t_future,
t_month = t_month_future),
interval = "confidence", level = 0.95)
# Combine actual data with forecast
forecast_2020_2021 <- data.frame(
Date = c(time(Arrivals_ts), time(future_dates)),
Actual = c(as.numeric(Arrivals_ts), rep(NA, length(future_dates))),
Forecast = c(rep(NA, length(Arrivals_ts)), predictions_future[,"fit"]),
Lower_CI = c(rep(NA, length(Arrivals_ts)), predictions_future[,"lwr"]),
Upper_CI = c(rep(NA, length(Arrivals_ts)), predictions_future[,"upr"]))
# Plot
ggplot(forecast_2020_2021, aes(x = Date)) +
geom_line(aes(y = Actual/1000000, color = "Actual")) +
geom_line(aes(y = Forecast/1000000, color = "Forecast")) +
geom_ribbon(aes(ymin = Lower_CI/1000000, ymax = Upper_CI/1000000), fill = "grey30", alpha = 0.4) +
labs(title = "Tourist Arrivals Forecast (2020-2021) with 95% CI",
y = "Tourist Arrivals (millions)") +
theme_minimal() +
scale_color_manual(values = c("Actual" = "darkblue", "Forecast" = "red"))
This linear model for monthly tourist exports data reveals a significant upward trend over time, with an estimated increase of 736.80 in exports per time unit (p < 0.05). The intercept, estimated at 13,864.95, represents the baseline level of exports. The high adjusted R^2 value of 0.9211 indicates that 92.11% of the variance in exports is explained by time, supporting a strong fit. Additionally, the significant F-statistic (2792, p < 0.05) affirms the model’s relevance. However, residuals display non-normality, lack independence, and are non-stationary, signaling potential limitations for forecasting or inference reliability in this model.
t <- time(Exports_ts)
t_centered <- t - mean(t)
# Linear
ex_lm <- lm(Exports_ts ~ t_centered)
summary(ex_lm)
Call:
lm(formula = Exports_ts ~ t_centered)
Residuals:
Min 1Q Median 3Q Max
-2551.6 -929.9 -169.9 701.6 3326.2
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 13864.95 80.50 172.23 <2e-16 ***
t_centered 736.80 13.94 52.84 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1247 on 238 degrees of freedom
Multiple R-squared: 0.9215, Adjusted R-squared: 0.9211
F-statistic: 2792 on 1 and 238 DF, p-value: < 2.2e-16ts.plot(cbind(Exports_ts/1000,ex_lm$fitted.values/1000),
col=c("darkblue", "darkred"), ylab = "Exports (billions $)",
main = "Original Data with Fitted Values - Linear", lwd =1.5)
legend("topleft", legend=c("Original Data", "Fitted Values"),
col=c("darkblue", "darkred"), lty=1:1,lwd=2)
checkresiduals(ex_lm$residuals)
Ljung-Box test
data: Residuals
Q* = 1464.4, df = 10, p-value < 2.2e-16
Model df: 0. Total lags used: 10check_stationarity(ex_lm$residuals)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| ADF Test | -3.005 | 0.15 | Non-Stationary |
| Phillips-Perron Test | -8.956 | 0.61 | Non-Stationary |
| KPSS Test | 0.515 | 0.04 | Non-Stationary |
| ERS Test | -0.496 | N/A | Non-Stationary |
check_normality_residuals(ex_lm)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| Shapiro-Wilk Test | 0.963 | 0 | Non-Normal |
| Anderson-Darling Test | 2.712 | 0 | Non-Normal |
| Jarque-Bera Test | 15.144 | 0.001 | Non-Normal |
| Lilliefors Test | 0.085 | 0 | Non-Normal |
| Skewness | 0.620 | N/A | Right-Skewed |
| Kurtosis | 3.010 | N/A | Mesokurtic |
In this quadratic model for monthly tourist exports, both the linear and quadratic time terms significantly predict exports, indicating a curved upward trend. The intercept, estimated at -5.899e+07, reflects the baseline export level when accounting for both terms. The quadratic term’s positive estimate (1.461 per time unit squared) suggests an accelerating increase over time. The high R^2 of 0.9305 implies that 93.05% of the variance in exports is explained by the model, and the substantial F-statistic (1602, p < 0.05) highlights its strong fit. Nonetheless, residuals remain non-normal, non-independent, and non-stationary, which may impact the model’s robustness for further interpretation or predictive use.
Note that the adjusted R^2 here is higher than that of the linear model alone, so we will keep this quadratic term in the following models
# Quadratic
ex_quad <- lm(Exports_ts ~ t_centered + I(t^2) )
summary(ex_quad)
Call:
lm(formula = Exports_ts ~ t_centered + I(t^2))
Residuals:
Min 1Q Median 3Q Max
-2504.8 -837.3 -230.1 968.2 2436.5
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.899e+07 1.024e+07 -5.762 2.57e-08 ***
t_centered -5.797e+04 1.019e+04 -5.691 3.71e-08 ***
I(t^2) 1.461e+01 2.534e+00 5.763 2.55e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1170 on 237 degrees of freedom
Multiple R-squared: 0.9311, Adjusted R-squared: 0.9305
F-statistic: 1602 on 2 and 237 DF, p-value: < 2.2e-16ts.plot(cbind(Exports_ts/1000,ex_quad$fitted.values/1000),
col=c("darkblue", "darkred"), ylab = "Exports (billions $)",
main = "Original Data with Fitted Values - Quadratic", lwd =1.5)
legend("topleft", legend=c("Original Data", "Fitted Values"),
col=c("darkblue", "darkred"), lty=1:1,lwd=2)
checkresiduals(ex_quad$residuals)
Ljung-Box test
data: Residuals
Q* = 1407.8, df = 10, p-value < 2.2e-16
Model df: 0. Total lags used: 10check_stationarity(ex_quad$residuals)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| ADF Test | -2.229 | 0.48 | Non-Stationary |
| Phillips-Perron Test | -7.588 | 0.68 | Non-Stationary |
| KPSS Test | 0.421 | 0.07 | Stationary |
| ERS Test | -0.611 | N/A | Non-Stationary |
check_normality_residuals(ex_quad)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| Shapiro-Wilk Test | 0.969 | 0 | Non-Normal |
| Anderson-Darling Test | 2.727 | 0 | Non-Normal |
| Jarque-Bera Test | 8.759 | 0.013 | Non-Normal |
| Lilliefors Test | 0.098 | 0 | Non-Normal |
| Skewness | 0.220 | N/A | Approximately Symmetric |
| Kurtosis | 2.180 | N/A | Platykurtic (Light Tails) |
This model incorporates both linear and quadratic time terms, along with monthly indicators, to capture trend and seasonality in the tourist exports data. Significant coefficients for the intercept, linear, and quadratic terms highlight a strong upward trend with a slight curvature. Monthly indicators are not statistically significant, suggesting that monthly seasonality has minimal impact on exports. The model’s adjusted R^2 of 0.9276 demonstrates strong explanatory power, supported by a high F-statistic (236.6, p < 0.05), indicating a good overall fit. However, residuals remain non-normal, non-independent, and non-stationary.
# Quadratic + Seasonal Trend
t_month <- factor(month(Masterdata$Time))
ex_month <- lm(Exports_ts ~ t_centered + I(t^2) + t_month)
summary(ex_month)
Call:
lm(formula = Exports_ts ~ t_centered + I(t^2) + t_month)
Residuals:
Min 1Q Median 3Q Max
-2483.1 -815.9 -244.8 970.2 2398.9
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.900e+07 1.045e+07 -5.645 4.92e-08 ***
t_centered -5.798e+04 1.040e+04 -5.575 7.01e-08 ***
I(t^2) 1.461e+01 2.587e+00 5.646 4.89e-08 ***
t_month2 1.255e+01 3.778e+02 0.033 0.974
t_month3 1.721e+01 3.778e+02 0.046 0.964
t_month4 -2.775e+01 3.779e+02 -0.073 0.942
t_month5 2.525e+01 3.779e+02 0.067 0.947
t_month6 7.344e+00 3.779e+02 0.019 0.985
t_month7 -7.287e+01 3.779e+02 -0.193 0.847
t_month8 -8.033e+01 3.779e+02 -0.213 0.832
t_month9 -1.009e+02 3.779e+02 -0.267 0.790
t_month10 -2.041e+02 3.780e+02 -0.540 0.590
t_month11 -2.157e+02 3.780e+02 -0.571 0.569
t_month12 -2.261e+02 3.780e+02 -0.598 0.550
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1195 on 226 degrees of freedom
Multiple R-squared: 0.9315, Adjusted R-squared: 0.9276
F-statistic: 236.6 on 13 and 226 DF, p-value: < 2.2e-16ts.plot(cbind(Exports_ts/1000,ex_month$fitted.values/1000 ),
col=c("darkblue", "darkred"), ylab = "Exports (billions $)",
main = "Original Data with Fitted Values - Quadratic + Seasonal", lwd =1.5)
legend("topleft", legend=c("Original Data", "Fitted Values"),
col=c("darkblue", "darkred"), lty=1:1,lwd=2)
checkresiduals(ex_month$residuals)
Ljung-Box test
data: Residuals
Q* = 1411, df = 10, p-value < 2.2e-16
Model df: 0. Total lags used: 10check_stationarity(ex_month$residuals)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| ADF Test | -2.284 | 0.46 | Non-Stationary |
| Phillips-Perron Test | -7.907 | 0.67 | Non-Stationary |
| KPSS Test | 0.420 | 0.07 | Stationary |
| ERS Test | -0.705 | N/A | Non-Stationary |
check_normality_residuals(ex_month)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| Shapiro-Wilk Test | 0.967 | 0 | Non-Normal |
| Anderson-Darling Test | 2.920 | 0 | Non-Normal |
| Jarque-Bera Test | 8.708 | 0.013 | Non-Normal |
| Lilliefors Test | 0.113 | 0 | Non-Normal |
| Skewness | 0.220 | N/A | Approximately Symmetric |
| Kurtosis | 2.180 | N/A | Platykurtic (Light Tails) |
This model includes linear and quadratic time components, as well as a cosine term to capture cyclical effects in the tourist exports data. The significant coefficients for the intercept, linear, and quadratic terms suggest a strong upward trend with a curvature, while the cosine term is not statistically significant, indicating that cyclical seasonality is not a prominent factor in this model. The model’s adjusted R^2 of 0.9304 indicates high explanatory power, supported by a substantial F-statistic (1066, p < 0.05), reflecting a good model fit. However, residuals remain non-normal, non-independent, and non-stationary.
# Linear + Cosine Trend
cosine_t <- cos(2 * pi * t / 12)
ex_cos <- lm(Exports_ts ~ t_centered + I(t^2) + cosine_t)
summary(ex_cos)
Call:
lm(formula = Exports_ts ~ t_centered + I(t^2) + cosine_t)
Residuals:
Min 1Q Median 3Q Max
-2418.1 -871.3 -246.4 984.6 2525.5
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.640e+07 1.080e+07 -5.223 3.86e-07 ***
t_centered -5.539e+04 1.074e+04 -5.156 5.34e-07 ***
I(t^2) 1.396e+01 2.673e+00 5.224 3.84e-07 ***
cosine_t 9.239e+01 1.213e+02 0.762 0.447
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1171 on 236 degrees of freedom
Multiple R-squared: 0.9313, Adjusted R-squared: 0.9304
F-statistic: 1066 on 3 and 236 DF, p-value: < 2.2e-16ts.plot(cbind(Exports_ts/1000,ex_cos$fitted.values/1000 ), col=c("darkblue", "darkred"),
ylab = "Exports (billions $)",
main = "Original Data with Fitted Values - Quadratic + Cosine", lwd =1.5)
legend("topleft", legend=c("Original Data", "Fitted Values"),
col=c("darkblue", "darkred"), lty=1:1,lwd=2)
checkresiduals(ex_cos$residuals)
Ljung-Box test
data: Residuals
Q* = 1404.9, df = 10, p-value < 2.2e-16
Model df: 0. Total lags used: 10check_stationarity(ex_cos$residuals)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| ADF Test | -2.327 | 0.44 | Non-Stationary |
| Phillips-Perron Test | -7.893 | 0.67 | Non-Stationary |
| KPSS Test | 0.421 | 0.07 | Stationary |
| ERS Test | -0.594 | N/A | Non-Stationary |
check_normality_residuals(ex_cos)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| Shapiro-Wilk Test | 0.969 | 0 | Non-Normal |
| Anderson-Darling Test | 2.727 | 0 | Non-Normal |
| Jarque-Bera Test | 9.600 | 0.008 | Non-Normal |
| Lilliefors Test | 0.096 | 0 | Non-Normal |
| Skewness | 0.260 | N/A | Approximately Symmetric |
| Kurtosis | 2.170 | N/A | Platykurtic (Light Tails) |
This model includes linear, quadratic, and cosine terms, along wih monthly indicators. Both the linear and quadratic terms are statistically significant (p < 0.05), suggesting that the trend follows a parabolic pattern. The cosine term and all monthly indicators, however, are not statistically significant (p > 0.05), indicating that it may not add meaningful cyclical variation. With an adjusted R^2of 0.9275, this model explains a substantial portion of the variance in the data. The F-statistic (219.3, p < 0.05) supports the model’s overall significance. However, residual issues persist, with non-normality, non-independence, and signs of non-stationarity.
# Quadratic + Seasonal + Cosine
ex_both <- lm(Exports_ts ~ t_centered + I(t^2) + cosine_t + t_month)
summary(ex_both)
Call:
lm(formula = Exports_ts ~ t_centered + I(t^2) + cosine_t + t_month)
Residuals:
Min 1Q Median 3Q Max
-2396.3 -851.8 -235.7 988.7 2487.8
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.641e+07 1.102e+07 -5.117 6.65e-07 ***
t_centered -5.541e+04 1.097e+04 -5.051 9.08e-07 ***
I(t^2) 1.397e+01 2.729e+00 5.118 6.61e-07 ***
cosine_t 9.229e+01 1.239e+02 0.745 0.457
t_month2 1.236e+01 3.782e+02 0.033 0.974
t_month3 1.685e+01 3.782e+02 0.045 0.965
t_month4 -2.823e+01 3.782e+02 -0.075 0.941
t_month5 2.468e+01 3.782e+02 0.065 0.948
t_month6 6.731e+00 3.783e+02 0.018 0.986
t_month7 -7.349e+01 3.783e+02 -0.194 0.846
t_month8 -8.092e+01 3.783e+02 -0.214 0.831
t_month9 -1.015e+02 3.783e+02 -0.268 0.789
t_month10 -2.045e+02 3.783e+02 -0.540 0.589
t_month11 -2.160e+02 3.784e+02 -0.571 0.569
t_month12 -2.262e+02 3.784e+02 -0.598 0.551
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1196 on 225 degrees of freedom
Multiple R-squared: 0.9317, Adjusted R-squared: 0.9275
F-statistic: 219.3 on 14 and 225 DF, p-value: < 2.2e-16ts.plot(cbind(Exports_ts/1000,ex_both$fitted.values/1000 ), col=c("darkblue", "darkred"),
ylab = "Exports (billions $)",
main = "Original Data with Fitted Values - Quadratic + Seasonal + Cosine", lwd =1.5)
legend("topleft", legend=c("Original Data", "Fitted Values"),
col=c("darkblue", "darkred"), lty=1:1,lwd=2)
checkresiduals(ex_both$residuals)
Ljung-Box test
data: Residuals
Q* = 1408, df = 10, p-value < 2.2e-16
Model df: 0. Total lags used: 10check_stationarity(ex_both$residuals)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| ADF Test | -2.381 | 0.42 | Non-Stationary |
| Phillips-Perron Test | -8.196 | 0.65 | Non-Stationary |
| KPSS Test | 0.420 | 0.07 | Stationary |
| ERS Test | -0.684 | N/A | Non-Stationary |
check_normality_residuals(ex_both)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| Shapiro-Wilk Test | 0.968 | 0 | Non-Normal |
| Anderson-Darling Test | 2.912 | 0 | Non-Normal |
| Jarque-Bera Test | 9.514 | 0.009 | Non-Normal |
| Lilliefors Test | 0.111 | 0 | Non-Normal |
| Skewness | 0.260 | N/A | Approximately Symmetric |
| Kurtosis | 2.170 | N/A | Platykurtic (Light Tails) |
Here, a quadratic model emerged as the best fit, with an R-squared of 0.93, an AIC of 4077.33, and a BIC of 4091.25. However, these metrics across all five models were relatively similar, suggesting no model captured the data particularly better than the others.
# Summary
ex_models <- list("Linear Model" = ex_lm,
"Quadratic Model" = ex_quad,
"Quadratic + Seasonal" = ex_month,
"Quadratic + Cosine" = ex_cos,
"Quadratic + Seasonal + Cosine" = ex_both)
model_summary_table(ex_models)| Model | Adjusted R² | AIC | BIC |
|---|---|---|---|
| Linear Model | 0.92 | 4106.8 | 4117.25 |
| Quadratic Model | 0.93 | 4077.33 | 4091.25 |
| Quadratic + Seasonal | 0.93 | 4097.83 | 4150.04 |
| Quadratic + Cosine | 0.93 | 4078.74 | 4096.14 |
| Quadratic + Seasonal + Cosine | 0.93 | 4099.24 | 4154.93 |
Additionally, the quadratic model’s residuals, like those of the arrivals model, failed to meet stationarity, normality, and independence assumptions, limiting the reliability of its forecast intervals. It also possesses problem of collinearity.
check_model(ex_both)
Here, we try to decompose the time series. As shown in the plots, while there seems to be signs of seasonality with frequency of 12 months, the magnitude of these seasonality signals are not so significant.
# Classical Seasonal Decomposition by Moving Averages
plot(decompose(Exports_ts/1000))
# STL Decomposition
plot(stl(Exports_ts/1000, s.window = 12))
Based on all of the poor statistics and test results above with this quadratic model, we know for certain that this model’s predictive power should not be deemed reliable and all predictions have to be interpreted with caution. Hence, we skipped the cross-validation test and proceeded with the prediction.
# Fit Model on Full Data and Forecast Two Years Ahead (with 95% CI)
ex_pred <- lm(Arrivals_ts ~ t + I(t^2))
future_dates <- ts(rep(NA, 24), start = c(2020, 1), frequency = 12)
t_future <- time(future_dates)
predictions_future <- predict(ex_pred,
newdata = data.frame(t = t_future,
`I(t^2)` = t_future^2),
interval = "confidence", level = 0.95)
forecast_2020_2021 <- data.frame(
Date = c(time(Arrivals_ts), time(future_dates)),
Actual = c(as.numeric(Arrivals_ts), rep(NA, length(future_dates))),
Forecast = c(rep(NA, length(Arrivals_ts)), predictions_future[,"fit"]),
Lower_CI = c(rep(NA, length(Arrivals_ts)), predictions_future[,"lwr"]),
Upper_CI = c(rep(NA, length(Arrivals_ts)), predictions_future[,"upr"]))
ggplot(forecast_2020_2021, aes(x = Date)) +
geom_line(aes(y = Actual / 1000000, color = "Actual")) +
geom_line(aes(y = Forecast / 1000000, color = "Forecast")) +
geom_ribbon(aes(ymin = Lower_CI / 1000000, ymax = Upper_CI / 1000000), fill = "grey30", alpha = 0.4) +
labs(title = "Tourist Arrivals Forecast (2020-2021) with 95% CI",
y = "Tourist Arrivals (millions)") +
theme_minimal() +
scale_color_manual(values = c("Actual" = "darkblue", "Forecast" = "red"))
This time series is definitely not stationary, based on the time series plot, the tests, and the ACF/PACF plots.
ggplot(Masterdata, aes(x = Time, y = Total)) +
geom_line(color = "darkblue", size = 1) +
ggtitle("US Monthly Tourist Arrivals (2000-2024)") +
xlab("Month") +
ylab("Tourist Arrivals") +
scale_x_date(date_labels = "%b %Y", date_breaks = "36 months")+
theme_light()+
scale_y_continuous(labels = scales::label_number(scale_cut = cut_short_scale()))
check_stationarity(Arrivals_ts)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| ADF Test | -5.456 | 0.01 | Stationary |
| Phillips-Perron Test | -113.997 | 0.01 | Stationary |
| KPSS Test | 4.286 | 0.01 | Non-Stationary |
| ERS Test | -1.360 | N/A | Non-Stationary |
two_plots(Arrivals_ts)
Auto.arima() suggests an ARIMA(3,0,2)(2,1,0)[12] with drift. From the outset, this model seems to fit the model well, with stationary and independent residuals. Yet, the residuals are not normal, suggesting poor reliability of forecasts. This step is only for reference, and we will attempt to build our only models.
# Auto.arima() - ARIMA(3,0,2)(2,1,0)[12] with drift
ar_initial <- auto.arima(Arrivals_ts)
ar_initialSeries: Arrivals_ts
ARIMA(3,0,2)(2,1,0)[12] with drift
Coefficients:
ar1 ar2 ar3 ma1 ma2 sar1 sar2 drift
0.8242 -0.2811 0.3483 -0.317 0.3689 -0.5835 -0.3680 13193.359
s.e. 0.2582 0.2985 0.1605 0.254 0.2163 0.0655 0.0682 6498.454
sigma^2 = 5.956e+10: log likelihood = -3151.28
AIC=6320.56 AICc=6321.38 BIC=6351.42ts.plot(cbind(Arrivals_ts/1000000,ar_initial$fitted/1000000),
col=c("darkblue", "darkred"), ylab = "Arrivals (millions)",
main = "Original Data with Fitted Values - Auto ARIMA", lwd =1.5)
legend("topleft", legend=c("Original Data", "Fitted Values"),
col=c("darkblue", "darkred"), lty=1:1,lwd=2)
checkresiduals(ar_initial)
Ljung-Box test
data: Residuals from ARIMA(3,0,2)(2,1,0)[12] with drift
Q* = 18.224, df = 17, p-value = 0.3748
Model df: 7. Total lags used: 24two_plots(ar_initial$residuals)
check_stationarity(ar_initial$residuals)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| ADF Test | -5.533 | 0.01 | Stationary |
| Phillips-Perron Test | -237.742 | 0.01 | Stationary |
| KPSS Test | 0.135 | 0.1 | Stationary |
| ERS Test | -7.112 | N/A | Stationary |
check_normality_residuals(ar_initial)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| Shapiro-Wilk Test | 0.967 | 0 | Non-Normal |
| Anderson-Darling Test | 2.200 | 0 | Non-Normal |
| Jarque-Bera Test | 51.371 | 0 | Non-Normal |
| Lilliefors Test | 0.074 | 0.003 | Non-Normal |
| Skewness | -0.230 | N/A | Approximately Symmetric |
| Kurtosis | 5.220 | N/A | Leptokurtic (Heavy Tails) |
Recall the ACF plot with a slow decay in both seasonal and normal lags, we will try both techniques here.
# Seasonal Differencing
sd_ar <- diff(Arrivals_ts,12)
diff_plots(sd_ar/1000000,"Seasonally Differenced Time Series")
check_stationarity(sd_ar)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| ADF Test | -3.026 | 0.14 | Non-Stationary |
| Phillips-Perron Test | -70.739 | 0.01 | Stationary |
| KPSS Test | 0.305 | 0.1 | Stationary |
| ERS Test | -4.092 | N/A | Stationary |
By only seasonally differencing, the series does not pass all stationarity tests.
# Normal Differencing
diff_ar <- diff(Arrivals_ts)
diff_plots(diff_ar, "Normally Differenced Time Series")
check_stationarity(diff_ar)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| ADF Test | -12.989 | 0.01 | Stationary |
| Phillips-Perron Test | -221.516 | 0.01 | Stationary |
| KPSS Test | 0.009 | 0.1 | Stationary |
| ERS Test | -5.811 | N/A | Stationary |
By only normally differencing, the series is already stationary, based on the “by-eye” test and othr statistical tests. Looking at the seasonal lags 12 and 24 decaying slowly in the ACF and the signicant lag 12 in the PACF, if we start fitting with only normal differencing, we will start with SAR 1. Note that the mean of residuals is zero, so we don’t need to include a drift.
# Normal & Seasonal Differencing
ar_nsd <- diff(diff(Arrivals_ts,12))
diff_plots(ar_nsd, "Normally & Seasonally Differenced Time Series")
check_stationarity(ar_nsd)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| ADF Test | -6.957 | 0.01 | Stationary |
| Phillips-Perron Test | -302.729 | 0.01 | Stationary |
| KPSS Test | 0.018 | 0.1 | Stationary |
| ERS Test | -2.911 | N/A | Stationary |
By differencing both ways, we achieved a stationary series, with smaller variance and less autocorrelations. Based on the ACF and PACF, a next move could be MA 1. Note that the mean of residuals is zero, so we don’t need to include a drift.
In this attempt, we will only focus on SARIMA models with \(d=1, D=0\).
First, we fit an ARIMA(0,1,0)(1,0,0)[12]. Based on the ACF and PACF, the best next move would be to add an AR1, or \(p=1\).
# ARIMA 010100
ar_010100 <- arima(Arrivals_ts, order=c(0,1,0),
seasonal=list(order=c(1,0,0), period=12))
ar_010100
Call:
arima(x = Arrivals_ts, order = c(0, 1, 0), seasonal = list(order = c(1, 0, 0),
period = 12))
Coefficients:
sar1
0.9001
s.e. 0.0238
sigma^2 estimated as 1.04e+11: log likelihood = -3380.49, aic = 6764.97checkresiduals(ar_010100)
Ljung-Box test
data: Residuals from ARIMA(0,1,0)(1,0,0)[12]
Q* = 109.87, df = 23, p-value = 2.663e-13
Model df: 1. Total lags used: 24two_plots(ar_010100$residuals)
check_stationarity(ar_010100$residuals)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| ADF Test | -7.658 | 0.01 | Stationary |
| Phillips-Perron Test | -310.831 | 0.01 | Stationary |
| KPSS Test | 0.015 | 0.1 | Stationary |
| ERS Test | -7.790 | N/A | Stationary |
check_normality_residuals(ar_010100)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| Shapiro-Wilk Test | 0.973 | 0 | Non-Normal |
| Anderson-Darling Test | 1.812 | 0 | Non-Normal |
| Jarque-Bera Test | 34.803 | 0 | Non-Normal |
| Lilliefors Test | 0.069 | 0.008 | Non-Normal |
| Skewness | -0.100 | N/A | Approximately Symmetric |
| Kurtosis | 4.850 | N/A | Leptokurtic (Heavy Tails) |
Here, we fit an ARIMA(0,1,1)(1,0,0)[12]. Based on the ACF and PACF, the best next move would be to add an SMA1, or \(Q=1\).
# ARIMA 011100
ar_011100 <- arima(Arrivals_ts, order=c(0,1,1),
seasonal=list(order=c(1,0,0), period=12))
ar_011100
Call:
arima(x = Arrivals_ts, order = c(0, 1, 1), seasonal = list(order = c(1, 0, 0),
period = 12))
Coefficients:
ma1 sar1
-0.4773 0.9175
s.e. 0.0585 0.0211
sigma^2 estimated as 8.338e+10: log likelihood = -3355.34, aic = 6716.69checkresiduals(ar_011100)
Ljung-Box test
data: Residuals from ARIMA(0,1,1)(1,0,0)[12]
Q* = 54.975, df = 22, p-value = 0.0001198
Model df: 2. Total lags used: 24two_plots(ar_011100$residuals)
check_stationarity(ar_011100$residuals)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| ADF Test | -6.288 | 0.01 | Stationary |
| Phillips-Perron Test | -251.548 | 0.01 | Stationary |
| KPSS Test | 0.022 | 0.1 | Stationary |
| ERS Test | -7.745 | N/A | Stationary |
check_normality_residuals(ar_011100)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| Shapiro-Wilk Test | 0.973 | 0 | Non-Normal |
| Anderson-Darling Test | 1.573 | 0 | Non-Normal |
| Jarque-Bera Test | 40.733 | 0 | Non-Normal |
| Lilliefors Test | 0.067 | 0.011 | Non-Normal |
| Skewness | 0.040 | N/A | Approximately Symmetric |
| Kurtosis | 5.020 | N/A | Leptokurtic (Heavy Tails) |
Here, we fit an ARIMA(0,1,1)(1,0,1)[12]. Based on the ACF, PACF, and other tests, we have finally found a model with independent and stationary residuals. All coefficients in the model are significant with small standard errors. However, the residuals are not normal; hence, we cannot trust the forecast intervals.
# ARIMA 011101
ar_011101 <- arima(Arrivals_ts, order=c(0,1,1),
seasonal=list(order=c(1,0,1), period=12), method="ML")
ar_011101
Call:
arima(x = Arrivals_ts, order = c(0, 1, 1), seasonal = list(order = c(1, 0, 1),
period = 12), method = "ML")
Coefficients:
ma1 sar1 sma1
-0.4569 0.9929 -0.6560
s.e. 0.0605 0.0036 0.0523
sigma^2 estimated as 5.924e+10: log likelihood = -3319.68, aic = 6647.36checkresiduals(ar_011101)
Ljung-Box test
data: Residuals from ARIMA(0,1,1)(1,0,1)[12]
Q* = 21.741, df = 21, p-value = 0.4146
Model df: 3. Total lags used: 24two_plots(ar_011101$residuals)
check_stationarity(ar_011101$residuals)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| ADF Test | -6.653 | 0.01 | Stationary |
| Phillips-Perron Test | -240.783 | 0.01 | Stationary |
| KPSS Test | 0.077 | 0.1 | Stationary |
| ERS Test | -7.521 | N/A | Stationary |
check_normality_residuals(ar_011101)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| Shapiro-Wilk Test | 0.982 | 0.004 | Non-Normal |
| Anderson-Darling Test | 1.064 | 0.008 | Non-Normal |
| Jarque-Bera Test | 20.178 | 0 | Non-Normal |
| Lilliefors Test | 0.052 | 0.115 | Normal |
| Skewness | -0.050 | N/A | Approximately Symmetric |
| Kurtosis | 4.420 | N/A | Leptokurtic (Heavy Tails) |
From the suggestions during differencing, we first fit an ARIMA(0,1,1)(0,1,0)[12]. Based on the ACF and PACF, a SMA 1 (\(Q=1\))should be recalled.
# ARIMA 011010
ar_011010 <- arima(Arrivals_ts, order=c(0,1,1), seasonal=list(order=c(0,1,0), period=12))
ar_011010
Call:
arima(x = Arrivals_ts, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 0),
period = 12))
Coefficients:
ma1
-0.4833
s.e. 0.0580
sigma^2 estimated as 8.707e+10: log likelihood = -3181.29, aic = 6366.59checkresiduals(ar_011010)
Ljung-Box test
data: Residuals from ARIMA(0,1,1)(0,1,0)[12]
Q* = 79.472, df = 23, p-value = 3.865e-08
Model df: 1. Total lags used: 24two_plots(ar_011010$residuals)
check_stationarity(ar_011010$residuals)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| ADF Test | -5.740 | 0.01 | Stationary |
| Phillips-Perron Test | -258.543 | 0.01 | Stationary |
| KPSS Test | 0.023 | 0.1 | Stationary |
| ERS Test | -7.880 | N/A | Stationary |
check_normality_residuals(ar_011010)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| Shapiro-Wilk Test | 0.963 | 0 | Non-Normal |
| Anderson-Darling Test | 2.302 | 0 | Non-Normal |
| Jarque-Bera Test | 66.311 | 0 | Non-Normal |
| Lilliefors Test | 0.072 | 0.004 | Non-Normal |
| Skewness | 0.030 | N/A | Approximately Symmetric |
| Kurtosis | 5.570 | N/A | Leptokurtic (Heavy Tails) |
Here, we fit an ARIMA(0,1,1)(0,1,1)[12]. Based on the ACF, PACF, and other statistical tests, we have found a model with independent and stationary residuals. However, the residuals are not normal; hence, we cannot trust the forecast intervals.
# ARIMA 011011
ar_011011 <- arima(Arrivals_ts, order=c(0,1,1),
seasonal=list(order=c(0,1,1), period=12))
ar_011011
Call:
arima(x = Arrivals_ts, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 1),
period = 12))
Coefficients:
ma1 sma1
-0.4600 -0.6602
s.e. 0.0609 0.0510
sigma^2 estimated as 5.947e+10: log likelihood = -3141.44, aic = 6288.87checkresiduals(ar_011011)
Ljung-Box test
data: Residuals from ARIMA(0,1,1)(0,1,1)[12]
Q* = 19.923, df = 22, p-value = 0.5879
Model df: 2. Total lags used: 24two_plots(ar_011011$residuals)
check_stationarity(ar_011011$residuals)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| ADF Test | -6.436 | 0.01 | Stationary |
| Phillips-Perron Test | -243.614 | 0.01 | Stationary |
| KPSS Test | 0.077 | 0.1 | Stationary |
| ERS Test | -7.635 | N/A | Stationary |
check_normality_residuals(ar_011011)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| Shapiro-Wilk Test | 0.981 | 0.002 | Non-Normal |
| Anderson-Darling Test | 1.374 | 0.001 | Non-Normal |
| Jarque-Bera Test | 19.836 | 0 | Non-Normal |
| Lilliefors Test | 0.073 | 0.003 | Non-Normal |
| Skewness | -0.010 | N/A | Approximately Symmetric |
| Kurtosis | 4.410 | N/A | Leptokurtic (Heavy Tails) |
Based on the summary table below, the best model appears to be ARIMA(0,1,1)(0,1,1)[12], with the lowest AIC and BIC values. Still, we recall that the residuals for this model are not normal, and thus we should not trust the forecast intervals with much reliability.
ar_sarima_models <- list("auto.arima()"= ar_initial,
"SARIMA(0,1,0)(1,0,0)[12] without drift" = ar_010100,
"SARIMA(0,1,1)(1,0,0)[12] without drift" = ar_011100,
"SARIMA(0,1,1)(1,0,1)[12] without drift" = ar_011101,
"SARIMA(0,1,1)(0,1,0)[12] without drift" = ar_011010,
"SARIMA(0,1,1)(0,1,1)[12] without drift" = ar_011011)
sarima_summary_table(ar_sarima_models)| Model | AIC | BIC |
|---|---|---|
| auto.arima() | 6320.56 | 6351.42 |
| SARIMA(0,1,0)(1,0,0)[12] without drift | 6764.97 | 6771.92 |
| SARIMA(0,1,1)(1,0,0)[12] without drift | 6716.69 | 6727.11 |
| SARIMA(0,1,1)(1,0,1)[12] without drift | 6647.36 | 6661.27 |
| SARIMA(0,1,1)(0,1,0)[12] without drift | 6366.59 | 6373.44 |
| SARIMA(0,1,1)(0,1,1)[12] without drift | 6288.87 | 6299.15 |
check_model(ar_011011)
After verifying that ARIMA(0,1,1)(0,1,1)[12] is also the best fit for the time period between 2000 and 2018, we used this model to predict the 2019 value and tested against the actual data. The result is astonishing with the fitted values following closely actual values. Moreover, an R^2 value of 0.9707 suggests that 97.07% of the variability in actual data can be explained by our model. We now want to leverage this model to predict the 2020 and 2021 hypothetical scenarios of no COVID.
#Setting up
train_data <- window(Arrivals_ts, end = c(2018, 12))
test_data <- window(Arrivals_ts, start = c(2019, 1), end = c(2019, 12))
forecast_2019 <- sarima.for(train_data, n.ahead = 12, 0, 1, 1, 0, 1, 1, 12, plot =FALSE)
predicted_2019 <- forecast_2019$pred
lower_ci_2019 <- forecast_2019$pred - 1.96 * forecast_2019$se
upper_ci_2019 <- forecast_2019$pred + 1.96 * forecast_2019$se
time_full_2019 <- data.frame(
Date = c(time(train_data), time(test_data)),
Actual = c(as.numeric(train_data), as.numeric(test_data)),
Predicted = c(rep(NA, length(train_data)), predicted_2019),
Lower_CI = c(rep(NA, length(train_data)), lower_ci_2019),
Upper_CI = c(rep(NA, length(train_data)), upper_ci_2019))
time_only_2019 <- data.frame(
Date = as.Date(time(test_data)),
Actual = test_data,
Predicted = predicted_2019,
Lower_CI = lower_ci_2019,
Upper_CI = upper_ci_2019)
# Plotting only 2019
ar_pred_2019_sarima <- ggplot(time_only_2019, aes(x = Date)) +
geom_line(aes(y = Actual / 1000000, color = "Actual")) +
geom_line(aes(y = Predicted / 1000000, color = "Predicted")) +
geom_ribbon(aes(ymin = Lower_CI / 1000000, ymax = Upper_CI / 1000000), fill = "grey30", alpha = 0.4) +
labs(title = "Actual and Predicted Arrivals (2019) with 95% CI",
y = "Tourist Arrivals (millions)") +
scale_x_date(date_labels = "%b %Y", date_breaks = "2 months") +
theme_minimal() +
scale_color_manual(values = c("Actual" = "darkblue", "Predicted" = "red"))
# Plotting full
ar_pred_full_sarima <- ggplot(time_full_2019, aes(x = Date)) +
geom_line(aes(y = Actual / 1000000, color = "Actual")) +
geom_line(aes(y = Predicted / 1000000, color = "Predicted")) +
geom_ribbon(aes(ymin = Lower_CI / 1000000, ymax = Upper_CI / 1000000), fill = "grey30", alpha = 0.4) +
labs(title = "Actual and Predicted Arrivals (2019) with 95% CI",
y = "Tourist Arrivals (millions)") +
theme_minimal() +
scale_color_manual(values = c("Actual" = "darkblue", "Predicted" = "red"))
grid.arrange(ar_pred_full_sarima, ar_pred_2019_sarima, nrow = 2)
# R^Squared
calculate_r2(as.numeric(predicted_2019), as.numeric(test_data))R-squared of the model is: 0.9707 Here, we can see that the model predicts monthly arrivals to be relatively stable, exhibiting historical seasonal trends with a slightly upward motion. This makes sense as without COVID-19, the US should have remained its attractiveness as an ideal place for travel, business, and education in 2020 and 2021.
# Run SARIMA forecast for 24 months ahead (2020 and 2021)
ar_sarima_for <- sarima.for(Arrivals_ts, n.ahead = 24, 0,1,1,0,1,1, S = 12, plot = FALSE)
forecasted_values <- ar_sarima_for$pred
lower_bound <- ar_sarima_for$pred - 1.96 * ar_sarima_for$se
upper_bound <- ar_sarima_for$pred + 1.96 * ar_sarima_for$se
time_forecast <- ts(rep(NA,24), start = c(2020, 1), end = c(2021, 12), frequency = 12)
forecast_ts <- ts(forecasted_values, start = c(2020, 1), frequency = 12)
lower_bound_ts <- ts(lower_bound, start = c(2020, 1), frequency = 12)
upper_bound_ts <- ts(upper_bound, start = c(2020, 1), frequency = 12)
combined_ts <- ts(c(Arrivals_ts, forecast_ts), start = start(Arrivals_ts), frequency = 12)
time_full_forecast <- data.frame(
Date = c(time(Arrivals_ts), time(time_forecast)),
Actual = c(as.numeric(Arrivals_ts), rep(NA, length(forecast_ts))),
Predicted = c(rep(NA, length(Arrivals_ts)), forecasted_values),
Lower_CI = c(rep(NA, length(Arrivals_ts)), lower_bound),
Upper_CI = c(rep(NA, length(Arrivals_ts)), upper_bound))
time_forecast_2020_2021 <- data.frame(
Date = as.Date(time(time_forecast)),
Predicted = forecasted_values,
Lower_CI = lower_bound,
Upper_CI = upper_bound)
# Plotting full dataset with forecasts
ar_pred_full_forecast <- ggplot(time_full_forecast, aes(x = Date)) +
geom_line(aes(y = Actual / 1000000, color = "Actual")) +
geom_line(aes(y = Predicted / 1000000, color = "Predicted")) +
geom_ribbon(aes(ymin = Lower_CI / 1000000, ymax = Upper_CI / 1000000), fill = "grey30", alpha = 0.4) +
labs(title = "Full Data with Forecasted Arrivals (2020-2021) with 95% CI",
y = "Tourist Arrivals (millions)") +
theme_minimal() +
scale_color_manual(values = c("Actual" = "darkblue", "Predicted" = "red"))
# Plotting only 2020 and 2021
ar_pred_2020_2021 <- ggplot(time_forecast_2020_2021, aes(x = Date)) +
geom_line(aes(y = Predicted / 1000000, color = "Predicted")) +
geom_ribbon(aes(ymin = Lower_CI / 1000000, ymax = Upper_CI / 1000000), fill = "grey30", alpha = 0.4) +
labs(title = "Forecasted Tourist Arrivals (2020-2021) with 95% CI",
y = "Tourist Arrivals (millions)") +
scale_x_date(date_labels = "%b %Y", date_breaks = "3 months") +
theme_minimal() +
scale_color_manual(values = c("Predicted" = "red"))
grid.arrange(ar_pred_full_forecast, ar_pred_2020_2021, nrow = 2)
This time series is definitely not stationary, based on the time series plot, the tests, and the ACF/PACF plots.
ggplot(Masterdata, aes(x = Time, y = Exports)) +
geom_line(color = "darkred", size = 1) +
ggtitle("US Monthly Tourist Exports, seasonally adjusted (2000-2024)") +
xlab("Month") +
ylab("Export Values (million $)") +
scale_x_date(date_labels = "%b %Y", date_breaks = "36 months")+
theme_light()+
scale_y_continuous(labels = scales::label_number(scale_cut = cut_short_scale()))
check_stationarity(Exports_ts)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| ADF Test | -3.005 | 0.15 | Non-Stationary |
| Phillips-Perron Test | -8.956 | 0.61 | Non-Stationary |
| KPSS Test | 4.746 | 0.01 | Non-Stationary |
| ERS Test | 1.206 | N/A | Non-Stationary |
two_plots(Exports_ts)
Auto.arima() suggests an ARIMA(2,1,2)(2,0,0)[12] with drift. From the outset, this model seems to fit the model well, with stationary and independent residuals. Yet, the residuals are not normal, suggesting poor reliability of forecasts. This step is only for reference, and we will attempt to build our only models.
# Auto.arima() - ARIMA(2,1,2)(2,0,0)[12] with drift
ex_initial <- auto.arima(Exports_ts)
ex_initialSeries: Exports_ts
ARIMA(2,1,2)(2,0,0)[12] with drift
Coefficients:
ar1 ar2 ma1 ma2 sar1 sar2 drift
-0.4818 -0.8840 0.5054 0.9604 -0.2794 -0.1773 47.1808
s.e. 0.0631 0.0615 0.0379 0.0378 0.0685 0.0769 12.5827
sigma^2 = 72552: log likelihood = -1674.01
AIC=3364.03 AICc=3364.65 BIC=3391.84ts.plot(cbind(Exports_ts/1000,ex_initial$fitted/1000),
col=c("darkblue", "darkred"), ylab = "Exports (billions $)",
main = "Original Data with Fitted Values - Auto ARIMA", lwd =1.5)
legend("topleft", legend=c("Original Data", "Fitted Values"),
col=c("darkblue", "darkred"), lty=1:1,lwd=2)
checkresiduals(ex_initial)
Ljung-Box test
data: Residuals from ARIMA(2,1,2)(2,0,0)[12] with drift
Q* = 17.561, df = 18, p-value = 0.4849
Model df: 6. Total lags used: 24two_plots(ex_initial$residuals)
check_stationarity(ex_initial$residuals)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| ADF Test | -5.506 | 0.01 | Stationary |
| Phillips-Perron Test | -226.347 | 0.01 | Stationary |
| KPSS Test | 0.297 | 0.1 | Stationary |
| ERS Test | -6.948 | N/A | Stationary |
check_normality_residuals(ex_initial)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| Shapiro-Wilk Test | 0.964 | 0 | Non-Normal |
| Anderson-Darling Test | 1.561 | 0.001 | Non-Normal |
| Jarque-Bera Test | 106.788 | 0 | Non-Normal |
| Lilliefors Test | 0.063 | 0.021 | Non-Normal |
| Skewness | -0.450 | N/A | Approximately Symmetric |
| Kurtosis | 6.140 | N/A | Leptokurtic (Heavy Tails) |
Recall the slow decay in the ACF plot, normal differencing is the first thing we should try. Luckily, it appears to be all that we need to achieve stationarity. Noting the non-zero mean of residuals, we will include a drift term in our models. Based on the ACF and PACF, the best next move appears to be a SMA 1 (\(Q=1\))
ex_diff <- diff(Exports_ts)
diff_plots(ex_diff,"Normally Differenced Time Series")
From the suggestions during differencing, we first fit an ARIMA(0,1,0)(0,0,1)[12]. Luckily, we immediately arrived at a good model with residuals showing stationarity and independence. However, the residuals are not normal; hence, we cannot trust the forecast intervals.
ex_010001 <- Arima(Exports_ts, order=c(0,1,0), seasonal=list(order=c(0,0,1), period=12), include.drift=TRUE)
ex_010001Series: Exports_ts
ARIMA(0,1,0)(0,0,1)[12] with drift
Coefficients:
sma1 drift
-0.2678 46.2844
s.e. 0.0725 13.1017
sigma^2 = 74133: log likelihood = -1678.6
AIC=3363.19 AICc=3363.29 BIC=3373.62checkresiduals(ex_010001)
Ljung-Box test
data: Residuals from ARIMA(0,1,0)(0,0,1)[12] with drift
Q* = 21.046, df = 23, p-value = 0.5783
Model df: 1. Total lags used: 24two_plots(ex_010001$residuals)
check_stationarity(ex_010001$residuals)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| ADF Test | -5.669 | 0.01 | Stationary |
| Phillips-Perron Test | -227.123 | 0.01 | Stationary |
| KPSS Test | 0.267 | 0.1 | Stationary |
| ERS Test | -6.824 | N/A | Stationary |
check_normality_residuals(ex_010001)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| Shapiro-Wilk Test | 0.960 | 0 | Non-Normal |
| Anderson-Darling Test | 1.577 | 0 | Non-Normal |
| Jarque-Bera Test | 115.892 | 0 | Non-Normal |
| Lilliefors Test | 0.083 | 0 | Non-Normal |
| Skewness | -0.430 | N/A | Approximately Symmetric |
| Kurtosis | 6.290 | N/A | Leptokurtic (Heavy Tails) |
From the ACF anf PACF plots during differencing, someone may suggest an SAR2 as well, i.e. (\(P=2\)). Hence, we fit an ARIMA(0,1,0)(2,0,0)[12]. Here, we also achieved stationarity and independence for residuals, though not normality. With two solid models, we wanted to see what would happen if we include more terms.
ex_010200 <- Arima(Exports_ts, order=c(0,1,0), seasonal=list(order=c(2,0,0), period=12), include.drift=TRUE)
ex_010200Series: Exports_ts
ARIMA(0,1,0)(2,0,0)[12] with drift
Coefficients:
sar1 sar2 drift
-0.2526 -0.1882 47.2346
s.e. 0.0665 0.0754 12.3734
sigma^2 = 73283: log likelihood = -1676.98
AIC=3361.95 AICc=3362.12 BIC=3375.86checkresiduals(ex_010200)
Ljung-Box test
data: Residuals from ARIMA(0,1,0)(2,0,0)[12] with drift
Q* = 22.37, df = 22, p-value = 0.438
Model df: 2. Total lags used: 24two_plots(ex_010200$residuals)
check_stationarity(ex_010200$residuals)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| ADF Test | -5.651 | 0.01 | Stationary |
| Phillips-Perron Test | -222.324 | 0.01 | Stationary |
| KPSS Test | 0.295 | 0.1 | Stationary |
| ERS Test | -6.752 | N/A | Stationary |
check_normality_residuals(ex_010200)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| Shapiro-Wilk Test | 0.964 | 0 | Non-Normal |
| Anderson-Darling Test | 1.557 | 0.001 | Non-Normal |
| Jarque-Bera Test | 103.015 | 0 | Non-Normal |
| Lilliefors Test | 0.077 | 0.001 | Non-Normal |
| Skewness | -0.410 | N/A | Approximately Symmetric |
| Kurtosis | 6.100 | N/A | Leptokurtic (Heavy Tails) |
ex_011200 <- Arima(Exports_ts, order=c(0,1,1), seasonal=list(order=c(2,0,0), period=12), include.drift=TRUE)
ex_011200Series: Exports_ts
ARIMA(0,1,1)(2,0,0)[12] with drift
Coefficients:
ma1 sar1 sar2 drift
0.0786 -0.2635 -0.1988 47.4203
s.e. 0.0602 0.0669 0.0752 13.1025
sigma^2 = 73033: log likelihood = -1676.13
AIC=3362.27 AICc=3362.52 BIC=3379.65checkresiduals(ex_011200)
Ljung-Box test
data: Residuals from ARIMA(0,1,1)(2,0,0)[12] with drift
Q* = 20.941, df = 21, p-value = 0.4626
Model df: 3. Total lags used: 24two_plots(ex_011200$residuals)
check_stationarity(ex_011200$residuals)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| ADF Test | -5.667 | 0.01 | Stationary |
| Phillips-Perron Test | -243.946 | 0.01 | Stationary |
| KPSS Test | 0.294 | 0.1 | Stationary |
| ERS Test | -6.871 | N/A | Stationary |
check_normality_residuals(ex_011200)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| Shapiro-Wilk Test | 0.966 | 0 | Non-Normal |
| Anderson-Darling Test | 1.527 | 0.001 | Non-Normal |
| Jarque-Bera Test | 84.386 | 0 | Non-Normal |
| Lilliefors Test | 0.070 | 0.007 | Non-Normal |
| Skewness | -0.320 | N/A | Approximately Symmetric |
| Kurtosis | 5.830 | N/A | Leptokurtic (Heavy Tails) |
ex_111200 <- Arima(Exports_ts, order=c(1,1,1), seasonal=list(order=c(2,0,0), period=12), include.drift=TRUE)
ex_111200Series: Exports_ts
ARIMA(1,1,1)(2,0,0)[12] with drift
Coefficients:
ar1 ma1 sar1 sar2 drift
0.4505 -0.3556 -0.2654 -0.1967 47.5588
s.e. 0.3912 0.4082 0.0669 0.0749 14.2048
sigma^2 = 73055: log likelihood = -1675.66
AIC=3363.32 AICc=3363.68 BIC=3384.18checkresiduals(ex_111200)
Ljung-Box test
data: Residuals from ARIMA(1,1,1)(2,0,0)[12] with drift
Q* = 20.388, df = 20, p-value = 0.4339
Model df: 4. Total lags used: 24two_plots(ex_111200$residuals)
check_stationarity(ex_111200$residuals)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| ADF Test | -5.863 | 0.01 | Stationary |
| Phillips-Perron Test | -240.199 | 0.01 | Stationary |
| KPSS Test | 0.272 | 0.1 | Stationary |
| ERS Test | -7.128 | N/A | Stationary |
check_normality_residuals(ex_111200)| Test | Test.Statistic | P.Value | Interpretation |
|---|---|---|---|
| Shapiro-Wilk Test | 0.964 | 0 | Non-Normal |
| Anderson-Darling Test | 1.663 | 0 | Non-Normal |
| Jarque-Bera Test | 91.445 | 0 | Non-Normal |
| Lilliefors Test | 0.066 | 0.013 | Non-Normal |
| Skewness | -0.290 | N/A | Approximately Symmetric |
| Kurtosis | 5.970 | N/A | Leptokurtic (Heavy Tails) |
Of all five models, our first two models appear to be the best fits to the data. However, while the second one (SARIMA(0,1,0)(2,0,0)[12] with drift) is more favorable in terms of AIC, the first one (SARIMA(0,1,0)(0,0,1)[12] with drift) is superior in terms of BIC. Now, as we know that the time series is seasonally adjusted, which means that seasonality should not play too much of a role in our prediction; hence, the notion that tourist exports from two years ago may still have an impact on that of this current month seems unplausible. Therefore, we believe that the first model, SARIMA(0,1,0)(0,0,1)[12] with drift, should be deemed as best fit. Again, as the residuals are not normal, we should not blindly trust our forecast intervals.
ex_sarima_models <- list("auto.arima() "= ex_initial,
"SARIMA(0,1,0)(0,0,1)[12] with drift" = ex_010001,
"SARIMA(0,1,0)(2,0,0)[12] with drift" = ex_010200,
"SARIMA(0,1,1)(2,0,0)[12] with drift" = ex_011200,
"SARIMA(1,1,1)(2,0,0)[12] with drift" = ex_111200)
sarima_summary_table(ex_sarima_models)| Model | AIC | BIC |
|---|---|---|
| auto.arima() | 3364.03 | 3391.84 |
| SARIMA(0,1,0)(0,0,1)[12] with drift | 3363.19 | 3373.62 |
| SARIMA(0,1,0)(2,0,0)[12] with drift | 3361.95 | 3375.86 |
| SARIMA(0,1,1)(2,0,0)[12] with drift | 3362.27 | 3379.65 |
| SARIMA(1,1,1)(2,0,0)[12] with drift | 3363.32 | 3384.18 |
check_model(ex_010001)
# Run SARIMA forecast for Exports_ts
ex_sarima_for <- sarima.for(Exports_ts, n.ahead = 24, 0,1,0,2,0,0, S = 12, plot = FALSE)
forecasted_values <- ex_sarima_for$pred
lower_bound <- ex_sarima_for$pred - 1.96 * ex_sarima_for$se
upper_bound <- ex_sarima_for$pred + 1.96 * ex_sarima_for$se
time_forecast <- time(Exports_ts)[length(Exports_ts)] + seq(1/12, 2, by=1/12)
ts_forecast <- ts(forecasted_values, start = c(2020, 1), frequency = 12)
combined_ts <- ts(c(Exports_ts, ts_forecast), start = start(Exports_ts), frequency = 12)
astsa::tsplot(combined_ts/1000, col = "darkblue", lwd = 1, type = "l",
ylab = "Exports (billions $)", main = "US Monthly Tourist Exports Forecast",
ylim = c(0, max(upper_bound/1000)))
lines(ts_forecast/1000, col = "red", lwd = 1.5) # Forecasted values
lines(ts(lower_bound, start = c(2020, 1), frequency = 12)/1000, lty = 3, col = "blue")
lines(ts(upper_bound, start = c(2020, 1), frequency = 12)/1000, lty = 3, col = "blue")
legend("topleft", legend = c("Observed", "Forecast", "95% Confidence Interval"),
col = c("darkblue", "red", "blue"), lty = c(1, 1, 3), lwd = c(1, 1.5, 1))