AUD/USD Exchange Rate Forecasting Analysis

MATH1318 Time Series Analysis - Final Project

Introduction

The Australian Dollar to United States Dollar (AUD/USD) exchange rate is one of the most actively traded currency pairs globally. This report presents a time series analysis of the AUD/USD monthly exchange rate from January 2010 to April 2026 using ARIMA and GARCH modelling.

Research Question: What will the AUD/USD exchange rate be for the next 10 months (May 2026 to February 2027)?

---

Section 0 - Load Packages

# Install and load all required packages
required_packages <- c("tidyverse", "lubridate", "forecast",
                       "tseries", "lmtest", "TSA",
                       "ggplot2", "gridExtra", "zoo", "moments")

for (pkg in required_packages) {
  if (!require(pkg, character.only = TRUE)) {
    install.packages(pkg)
    library(pkg, character.only = TRUE)
  }
}

---

Section 1 - Load and Prepare Data

# Set working directory — update this path to where your f11-data.csv is saved
setwd("C:/Users/Aleena/Downloads/")

# Load the CSV — skip the 10 metadata rows at the top of the RBA file
df_raw <- read.csv("f11-data.csv",
                   skip             = 10,
                   header           = TRUE,
                   stringsAsFactors = FALSE)

# Rename the first two columns to Date and AUD_USD
colnames(df_raw)[1] <- "Date"
colnames(df_raw)[2] <- "AUD_USD"

# Parse dates, coerce to numeric, remove NAs, keep only Date and AUD_USD columns
df <- df_raw %>%
  mutate(
    Date    = as.Date(Date, format = "%d-%b-%y"),
    AUD_USD = as.numeric(AUD_USD)
  ) %>%
  filter(!is.na(Date), !is.na(AUD_USD)) %>%
  select(Date, AUD_USD)

# Verify data loaded correctly
cat("=== Data loaded successfully ===\n")

## === Data loaded successfully ===

cat("Observations :", nrow(df), "\n")

## Observations : 196

cat("Date range   :", format(min(df$Date), "%b %Y"),
    "to", format(max(df$Date), "%b %Y"), "\n")

## Date range   : Jan 2010 to Apr 2026

cat("Missing values:", sum(is.na(df$AUD_USD)), "\n")

## Missing values: 0

# Create monthly time series object starting January 2010
ts_audusd <- ts(df$AUD_USD, start = c(2010, 1), frequency = 12)

---

Section 2 - Descriptive Analysis

Summary Statistics

# Descriptive statistics for the AUD/USD series
cat("=== Summary Statistics - AUD/USD ===\n")

## === Summary Statistics - AUD/USD ===

summary(ts_audusd)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.6175  0.6893  0.7462  0.7897  0.8947  1.0954

cat("Standard deviation:", round(sd(ts_audusd), 4), "\n")

## Standard deviation: 0.1337

cat("Variance          :", round(var(ts_audusd), 4), "\n")

## Variance          : 0.0179

cat("Skewness          :", round(moments::skewness(ts_audusd), 4), "\n")

## Skewness          : 0.8422

Time Series Plot

# Time series plot reveals trend, COVID-19 shock (2020), and volatility clustering
autoplot(ts_audusd) +
  geom_line(colour = "#2563EB", linewidth = 0.8) +
  labs(
    title    = "Figure 1. AUD/USD Monthly Exchange Rate (Jan 2010 to Apr 2026)",
    subtitle = "Source: Reserve Bank of Australia - F11 Historical Exchange Rates",
    x        = "Year",
    y        = "AUD/USD Rate"
  ) +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face = "bold"))

Figure 1. AUD/USD Monthly Exchange Rate (Jan 2010 to Apr 2026). Source: Reserve Bank of Australia F11.

Classical Decomposition

# Additive decomposition separates trend, seasonal, and remainder components
decomp <- decompose(ts_audusd, type = "additive")
autoplot(decomp) +
  labs(title = "Figure 2. Classical Additive Decomposition - AUD/USD") +
  theme_minimal()

Figure 2. Classical additive decomposition of the AUD/USD series.

Seasonal Subseries Plot

# Flat monthly means indicate no systematic seasonal pattern
monthplot(ts_audusd,
          main = "Figure 3. Seasonal Subseries Plot - AUD/USD",
          ylab = "AUD/USD",
          xlab = "Month")

Figure 3. Seasonal subseries plot showing monthly means across years.

ACF and PACF of Original Series

# Slow decay in ACF is characteristic of a non-stationary series with a unit root
par(mfrow = c(1, 2))
acf(ts_audusd,
    lag.max = 48,
    main    = "Figure 4a. ACF - Original Series")
pacf(ts_audusd,
     lag.max = 48,
     main    = "Figure 4b. PACF - Original Series")

Figure 4. ACF and PACF of the original AUD/USD series. Slow ACF decay confirms non-stationarity.

par(mfrow = c(1, 1))

---

Section 3 - Stationarity Testing

ADF and KPSS Tests on Original Series

# ADF test H0: series is non-stationary. p > 0.05 = fail to reject = non-stationary
cat("=== ADF Test - Original Series ===\n")

## === ADF Test - Original Series ===

adf_orig <- adf.test(ts_audusd)
print(adf_orig)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts_audusd
## Dickey-Fuller = -2.1039, Lag order = 5, p-value = 0.5326
## alternative hypothesis: stationary

# KPSS test H0: series is stationary. p < 0.05 = reject H0 = non-stationary
cat("\n=== KPSS Test - Original Series ===\n")

## 
## === KPSS Test - Original Series ===

kpss_orig <- kpss.test(ts_audusd)
print(kpss_orig)

## 
##  KPSS Test for Level Stationarity
## 
## data:  ts_audusd
## KPSS Level = 3.1395, Truncation lag parameter = 4, p-value = 0.01

First Differencing and Re-testing

# Apply first difference to remove the non-stationarity identified above
ts_diff <- diff(ts_audusd, differences = 1)

# ADF on differenced series: p < 0.05 confirms stationarity achieved
cat("=== ADF Test - First Differenced Series ===\n")

## === ADF Test - First Differenced Series ===

adf_diff <- adf.test(ts_diff)
print(adf_diff)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts_diff
## Dickey-Fuller = -4.9418, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

# KPSS on differenced series: p > 0.05 confirms stationarity
cat("\n=== KPSS Test - First Differenced Series ===\n")

## 
## === KPSS Test - First Differenced Series ===

kpss_diff <- kpss.test(ts_diff)
print(kpss_diff)

## 
##  KPSS Test for Level Stationarity
## 
## data:  ts_diff
## KPSS Level = 0.095709, Truncation lag parameter = 4, p-value = 0.1

Differenced Series Plot

# Differenced series should have no trend and constant variance
autoplot(ts_diff) +
  geom_line(colour = "#0D9488", linewidth = 0.7) +
  labs(title    = "Figure 5. AUD/USD - First Differenced Series",
       subtitle = "Stationary: mean near zero, no visible trend",
       x        = "Year",
       y        = "Delta AUD/USD") +
  theme_minimal()

Figure 5. First differenced AUD/USD series. Fluctuation around zero confirms stationarity (d=1).

ACF and PACF of Differenced Series

# ACF single spike at lag 1 then cuts off = MA(1) = q=1
# PACF no significant spikes = no AR component = p=0
# Together these identify ARIMA(0,1,1) as the candidate model
par(mfrow = c(1, 2))
acf(ts_diff,
    lag.max = 36,
    main    = "Figure 6a. ACF - Differenced Series (identifies q)")
pacf(ts_diff,
     lag.max = 36,
     main    = "Figure 6b. PACF - Differenced Series (identifies p)")

Figure 6. ACF and PACF of the first differenced series. ACF cuts off at lag 1 suggesting q=1; PACF shows no spikes suggesting p=0.

par(mfrow = c(1, 1))

---

Section 4 - Seasonality Assessment

ACF at Seasonal Lags

# No spikes at lags 12, 24, 36 = no seasonal autocorrelation = SARIMA not needed
acf(as.vector(ts_audusd),
    lag.max = 48,
    main    = "Figure 7. ACF to Lag 48 - Checking for Seasonal Spikes at 12, 24, 36")

Figure 7. ACF to lag 48. No significant spikes at lags 12, 24, or 36 indicate no seasonal component.

Seasonal Subseries Plot

# Monthly means that are approximately equal across all 12 months = no seasonality
monthplot(ts_audusd,
          main = "Figure 8. Seasonal Subseries Plot - AUD/USD",
          ylab = "AUD/USD",
          xlab = "Month")

Figure 8. Seasonal subseries plot. Flat monthly means confirm absence of seasonal pattern.

Residuals Approach - Step 1

# Fit seasonal-difference-only model and inspect residuals for seasonal pattern
m_seas <- arima(ts_audusd,
                order    = c(0, 0, 0),
                seasonal = list(order = c(0, 1, 0), period = 12))
res_seas <- residuals(m_seas)

# Time plot of residuals after seasonal difference
plot(res_seas,
     main = "Figure 9. Time Series Plot of Residuals after Seasonal Difference",
     ylab = "Residuals",
     xlab = "Time")
abline(h = 0, col = "red", lty = 2)

Figure 9. Residuals after seasonal differencing. No seasonal structure visible.

# No significant spikes at seasonal lags confirms no seasonal component
par(mfrow = c(1, 2))
acf(res_seas,
    lag.max = 36,
    main    = "Figure 10a. ACF of Residuals after Seasonal Difference")
pacf(res_seas,
     lag.max = 36,
     main    = "Figure 10b. PACF of Residuals after Seasonal Difference")

Figure 10. ACF and PACF of residuals after seasonal differencing. No seasonal autocorrelation remains.

par(mfrow = c(1, 1))

Residuals Approach - Step 2

# Add ordinary difference to remove any remaining trend in residuals
m_seas_diff <- arima(ts_audusd,
                     order    = c(0, 1, 0),
                     seasonal = list(order = c(0, 1, 0), period = 12))
res_seas_diff <- residuals(m_seas_diff)

par(mfrow = c(1, 2))
acf(res_seas_diff,
    lag.max = 36,
    main    = "Figure 11a. ACF after Ordinary and Seasonal Difference")
pacf(res_seas_diff,
     lag.max = 36,
     main    = "Figure 11b. PACF after Ordinary and Seasonal Difference")

Figure 11. ACF and PACF of residuals after ordinary and seasonal differencing.

par(mfrow = c(1, 1))

Seasonality Conclusion

# Summary of seasonality assessment
cat("=== Seasonality Assessment ===\n")

## === Seasonality Assessment ===

cat("Tools used: ACF to lag 48, seasonal subseries plot, residuals approach\n\n")

## Tools used: ACF to lag 48, seasonal subseries plot, residuals approach

cat("Findings:\n")

## Findings:

cat("  - ACF shows no significant spikes at lags 12, 24, 36\n")

##   - ACF shows no significant spikes at lags 12, 24, 36

cat("  - Seasonal subseries plot shows approximately flat monthly means\n")

##   - Seasonal subseries plot shows approximately flat monthly means

cat("  - Residual ACF/PACF after seasonal difference shows no seasonal pattern\n\n")

##   - Residual ACF/PACF after seasonal difference shows no seasonal pattern

cat("Conclusion:\n")

## Conclusion:

cat("  No seasonal component detected in the AUD/USD exchange rate series\n")

##   No seasonal component detected in the AUD/USD exchange rate series

cat("  SARIMA modelling is NOT appropriate for this series\n")

##   SARIMA modelling is NOT appropriate for this series

cat("  Proceeding with ARIMA (mean model) + GARCH (variance model)\n")

##   Proceeding with ARIMA (mean model) + GARCH (variance model)

---

Section 5 - ARIMA Model Specification and Fitting

auto.arima Benchmark

# auto.arima used as a benchmark to confirm the ACF/PACF-based identification
cat("=== auto.arima() Benchmark ===\n")

## === auto.arima() Benchmark ===

auto_model <- auto.arima(ts_audusd,
                          seasonal      = FALSE,
                          stepwise      = FALSE,
                          approximation = FALSE)
summary(auto_model)

## Series: ts_audusd 
## ARIMA(0,1,1) 
## 
## Coefficients:
##           ma1
##       -0.1363
## s.e.   0.0740
## 
## sigma^2 = 0.0006457:  log likelihood = 439.95
## AIC=-875.9   AICc=-875.83   BIC=-869.35
## 
## Training set error measures:
##                        ME       RMSE        MAE        MPE    MAPE      MASE
## Training set -0.001075096 0.02528141 0.01935988 -0.1881601 2.43201 0.3296916
##                     ACF1
## Training set 0.002059628

Candidate Model Comparison

# Fit six candidate ARIMA models covering a range of p and q orders
m1 <- Arima(ts_audusd, order = c(1, 1, 1))
m2 <- Arima(ts_audusd, order = c(1, 1, 0))
m3 <- Arima(ts_audusd, order = c(0, 1, 1))   # identified from ACF/PACF
m4 <- Arima(ts_audusd, order = c(2, 1, 1))
m5 <- Arima(ts_audusd, order = c(1, 1, 2))
m6 <- Arima(ts_audusd, order = c(2, 1, 2))

# Compare all models by AIC, BIC, and log-likelihood — lower AIC/BIC is better
model_comparison <- data.frame(
  Model  = c("ARIMA(1,1,1)", "ARIMA(1,1,0)", "ARIMA(0,1,1)",
             "ARIMA(2,1,1)", "ARIMA(1,1,2)", "ARIMA(2,1,2)"),
  AIC    = c(AIC(m1), AIC(m2), AIC(m3), AIC(m4), AIC(m5), AIC(m6)),
  BIC    = c(BIC(m1), BIC(m2), BIC(m3), BIC(m4), BIC(m5), BIC(m6)),
  LogLik = c(as.numeric(logLik(m1)), as.numeric(logLik(m2)),
             as.numeric(logLik(m3)), as.numeric(logLik(m4)),
             as.numeric(logLik(m5)), as.numeric(logLik(m6)))
) %>%
  mutate(across(where(is.numeric), ~ round(.x, 3))) %>%
  arrange(AIC)

cat("=== Model Comparison Table (sorted by AIC) ===\n")

## === Model Comparison Table (sorted by AIC) ===

print(model_comparison)

##          Model      AIC      BIC  LogLik
## 1 ARIMA(0,1,1) -875.897 -869.351 439.948
## 2 ARIMA(1,1,0) -875.637 -869.091 439.819
## 3 ARIMA(1,1,1) -874.304 -864.485 440.152
## 4 ARIMA(2,1,1) -872.310 -859.218 440.155
## 5 ARIMA(1,1,2) -872.309 -859.217 440.155
## 6 ARIMA(2,1,2) -871.539 -855.174 440.770

cat("\nBest model by AIC:", model_comparison$Model[1], "\n")

## 
## Best model by AIC: ARIMA(0,1,1)

Best Model Parameter Estimates

# Select ARIMA(0,1,1) — lowest AIC and BIC, confirmed by ACF/PACF analysis
best_arima <- m3

cat("=== Best ARIMA Model Summary ===\n")

## === Best ARIMA Model Summary ===

summary(best_arima)

## Series: ts_audusd 
## ARIMA(0,1,1) 
## 
## Coefficients:
##           ma1
##       -0.1363
## s.e.   0.0740
## 
## sigma^2 = 0.0006457:  log likelihood = 439.95
## AIC=-875.9   AICc=-875.83   BIC=-869.35
## 
## Training set error measures:
##                        ME       RMSE        MAE        MPE    MAPE      MASE
## Training set -0.001075096 0.02528141 0.01935988 -0.1881601 2.43201 0.3296916
##                     ACF1
## Training set 0.002059628

# CSS estimation method (as taught in Module 6)
cat("\n=== CSS Parameter Estimates ===\n")

## 
## === CSS Parameter Estimates ===

best_css <- arima(ts_audusd, order = c(0, 1, 1), method = "CSS")
coeftest(best_css)

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)  
## ma1 -0.137043   0.074149 -1.8482  0.06457 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# ML estimation method for comparison (as taught in Module 6)
cat("\n=== ML Parameter Estimates ===\n")

## 
## === ML Parameter Estimates ===

best_ml <- arima(ts_audusd, order = c(0, 1, 1), method = "ML")
coeftest(best_ml)

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)  
## ma1 -0.136298   0.073967 -1.8427  0.06537 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

---

Section 6 - ARIMA Diagnostic Checking

Standardised Residuals Plot

# Standardised residuals should be randomly distributed around zero
plot(rstandard(best_arima),
     ylab = "Standardised Residuals",
     type = "o",
     main = "Figure 12. Standardised Residuals - ARIMA(0,1,1)")
abline(h = 0)

Figure 12. Standardised residuals from ARIMA(0,1,1). Random fluctuation around zero is consistent with white noise assumption.

ACF of Standardised Residuals

# All ACF values within confidence bounds = residuals are uncorrelated = white noise
acf(as.vector(rstandard(best_arima)),
    lag.max = 36,
    main    = "Figure 13. ACF of Standardised Residuals - ARIMA(0,1,1)")

Figure 13. ACF of standardised residuals. All values within 95% bounds confirm no remaining autocorrelation.

Histogram and Q-Q Plot

# Histogram and QQ-plot to assess normality of residuals
par(mfrow = c(1, 2))
hist(rstandard(best_arima),
     xlab = "Standardised Residuals",
     main = "Figure 14a. Histogram of Residuals")
qqnorm(rstandard(best_arima),
       main = "Figure 14b. Q-Q Plot of Residuals")
qqline(rstandard(best_arima), col = "red", lwd = 2)

Figure 14. Histogram and Q-Q plot of standardised residuals showing approximate normality with slight fat tails.

par(mfrow = c(1, 1))

Ljung-Box and Shapiro-Wilk Tests

# Ljung-Box test H0: residuals are white noise. p > 0.05 = PASS
cat("=== Ljung-Box Test on Standardised Residuals ===\n")

## === Ljung-Box Test on Standardised Residuals ===

lb_test <- Box.test(rstandard(best_arima),
                    lag   = 20,
                    type  = "Ljung-Box",
                    fitdf = 0)
print(lb_test)

## 
##  Box-Ljung test
## 
## data:  rstandard(best_arima)
## X-squared = 22.744, df = 20, p-value = 0.3015

# Shapiro-Wilk test H0: residuals are normally distributed. p < 0.05 = slight fat tails
cat("\n=== Shapiro-Wilk Normality Test ===\n")

## 
## === Shapiro-Wilk Normality Test ===

shapiro.test(rstandard(best_arima))

## 
##  Shapiro-Wilk normality test
## 
## data:  rstandard(best_arima)
## W = 0.98434, p-value = 0.02807

Overfitting Check

# Fit ARIMA(0,1,2) to confirm ARIMA(0,1,1) is not under-specified
# If MA(2) is insignificant and AIC is worse, the simpler model is confirmed
cat("=== Overfitting Check - ARIMA(0,1,2) ===\n")

## === Overfitting Check - ARIMA(0,1,2) ===

m_overfit <- arima(ts_audusd, order = c(0, 1, 2), method = "ML")
coeftest(m_overfit)

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)  
## ma1 -0.133703   0.071533 -1.8691  0.06161 .
## ma2 -0.039283   0.074727 -0.5257  0.59910  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

AIC(best_arima, m_overfit)

---

Section 7 - ARCH Effects Testing

ACF and PACF of Squared and Absolute Residuals

# Significant autocorrelation in squared/absolute residuals = volatility clustering
library(TSA)
arima_resid <- residuals(best_arima)

par(mfrow = c(2, 2))
acf(arima_resid^2,     lag.max = 36, main = "Figure 15a. ACF - Squared Residuals")
pacf(arima_resid^2,    lag.max = 36, main = "Figure 15b. PACF - Squared Residuals")
acf(abs(arima_resid),  lag.max = 36, main = "Figure 15c. ACF - Absolute Residuals")
pacf(abs(arima_resid), lag.max = 36, main = "Figure 15d. PACF - Absolute Residuals")

Figure 15. ACF and PACF of squared and absolute ARIMA residuals. Significant spikes indicate volatility clustering (ARCH effects).

par(mfrow = c(1, 1))

McLeod-Li Test

# McLeod-Li test is the formal test for ARCH effects taught in Module 9
# p < 0.05 at most lags = significant ARCH effects = GARCH is justified
McLeod.Li.test(y    = arima_resid,
               main = "Figure 16. P-values of McLeod-Li Test - AUD/USD ARIMA Residuals")

Figure 16. P-values of McLeod-Li test. Values below the 5% line confirm ARCH effects are present, justifying GARCH modelling.

# Formal conclusion
cat("\n=== McLeod-Li Test Conclusion ===\n")

## 
## === McLeod-Li Test Conclusion ===

ml_test <- McLeod.Li.test(y = arima_resid, plot = FALSE)
if (any(ml_test$p.values < 0.05)) {
  cat("ARCH effects detected -> GARCH modelling is justified\n")
} else {
  cat("No significant ARCH effects detected\n")
}

## ARCH effects detected -> GARCH modelling is justified

---

Section 8 - GARCH Modelling

Log Returns Plot

# Compute percentage log returns for GARCH modelling as per Module 9
r_audusd <- diff(log(ts_audusd)) * 100

# Log returns plot should show volatility clustering
plot(r_audusd, type = "l",
     main = "Figure 17. AUD/USD Monthly Log Returns (%)",
     ylab = "Log Return (%)",
     xlab = "Time")
abline(h = 0)

Figure 17. AUD/USD monthly log returns. Volatility clustering visible with alternating calm and turbulent periods.

ACF/PACF for GARCH Order Identification

# ARMA pattern in squared returns identifies p and max(p,q) for GARCH model
par(mfrow = c(2, 2))
acf(r_audusd^2,      lag.max = 36, main = "Figure 18a. ACF - Squared Returns")
pacf(r_audusd^2,     lag.max = 36, main = "Figure 18b. PACF - Squared Returns")
acf(abs(r_audusd),   lag.max = 36, main = "Figure 18c. ACF - Absolute Returns")
pacf(abs(r_audusd),  lag.max = 36, main = "Figure 18d. PACF - Absolute Returns")

Figure 18. ACF and PACF of squared and absolute returns used to identify GARCH orders.

par(mfrow = c(1, 1))

GARCH Model Fitting and Comparison

# Fit GARCH(1,1) using tseries::garch — exact function used in Module 9 notes
m_garch11 <- garch(r_audusd, order = c(1, 1), trace = FALSE)
cat("=== GARCH(1,1) Model ===\n")

## === GARCH(1,1) Model ===

summary(m_garch11)

## 
## Call:
## garch(x = r_audusd, order = c(1, 1), trace = FALSE)
## 
## Model:
## GARCH(1,1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.18899 -0.66815 -0.04892  0.62646  2.56631 
## 
## Coefficient(s):
##     Estimate  Std. Error  t value Pr(>|t|)    
## a0   0.34822     0.65158    0.534    0.593    
## a1   0.02556     0.03199    0.799    0.424    
## b1   0.93356     0.08830   10.573   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Diagnostic Tests:
##  Jarque Bera Test
## 
## data:  Residuals
## X-squared = 0.62537, df = 2, p-value = 0.7315
## 
## 
##  Box-Ljung test
## 
## data:  Squared.Residuals
## X-squared = 0.21347, df = 1, p-value = 0.6441

# Fit GARCH(2,2) as comparison model
m_garch22 <- garch(r_audusd, order = c(2, 2), trace = FALSE)
cat("\n=== GARCH(2,2) Model ===\n")

## 
## === GARCH(2,2) Model ===

summary(m_garch22)

## 
## Call:
## garch(x = r_audusd, order = c(2, 2), trace = FALSE)
## 
## Model:
## GARCH(2,2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.03062 -0.64442 -0.05196  0.60150  2.30344 
## 
## Coefficient(s):
##     Estimate  Std. Error  t value Pr(>|t|)
## a0 7.611e+00   5.545e+00    1.373    0.170
## a1 1.246e-02   5.377e-02    0.232    0.817
## a2 1.580e-01   1.856e-01    0.852    0.394
## b1 3.004e-13   4.440e-01    0.000    1.000
## b2 3.684e-02   5.530e-01    0.067    0.947
## 
## Diagnostic Tests:
##  Jarque Bera Test
## 
## data:  Residuals
## X-squared = 0.7716, df = 2, p-value = 0.6799
## 
## 
##  Box-Ljung test
## 
## data:  Squared.Residuals
## X-squared = 0.077206, df = 1, p-value = 0.7811

# AIC comparison — GARCH(2,2) has lower AIC but all its parameters are insignificant
cat("\n=== GARCH Model AIC Comparison ===\n")

## 
## === GARCH Model AIC Comparison ===

cat("GARCH(1,1) AIC:", round(AIC(m_garch11), 4), "\n")

## GARCH(1,1) AIC: 991.7976

cat("GARCH(2,2) AIC:", round(AIC(m_garch22), 4), "\n")

## GARCH(2,2) AIC: 988.9276

cat("GARCH(1,1) selected: lower AIC is outweighed by all GARCH(2,2) params being insignificant\n")

## GARCH(1,1) selected: lower AIC is outweighed by all GARCH(2,2) params being insignificant

# Store best GARCH model
best_garch <- m_garch11

# Validity check: alpha1 + beta1 must be < 1 for mean-reverting volatility
garch_coef <- coef(best_garch)
cat("\nalpha1 + beta1 =",
    round(garch_coef["a1"] + garch_coef["b1"], 4),
    "-- must be < 1 for mean-reverting volatility\n")

## 
## alpha1 + beta1 = 0.9591 -- must be < 1 for mean-reverting volatility

GARCH Diagnostic Plots

# Standardised residuals from GARCH model should appear random
plot(residuals(best_garch), type = "h",
     ylab = "Standardised Residuals",
     main = "Figure 19. Standardised Residuals from GARCH(1,1)")

Figure 19. Standardised residuals from GARCH(1,1). Random pattern with no structure confirms model adequacy.

# Q-Q plot of GARCH standardised residuals
qqnorm(residuals(best_garch),
       main = "Figure 20. Q-Q Plot - Standardised Residuals from GARCH(1,1)")
qqline(residuals(best_garch))

Figure 20. Q-Q plot of GARCH standardised residuals. Close alignment with diagonal confirms approximate normality.

# No significant autocorrelation in squared residuals = ARCH effects removed
acf(residuals(best_garch)^2,
    na.action = na.omit,
    main = "Figure 21. ACF of Squared Standardised Residuals - GARCH(1,1)")

Figure 21. ACF of squared standardised GARCH residuals. No significant spikes confirm no remaining ARCH effects.

# Manual p-value plot equivalent to LBQPlot from Module 9
res_garch <- na.omit(residuals(best_garch))
lags      <- 1:24
pvalues   <- sapply(lags, function(lag) {
  Box.test(res_garch^2, lag = lag, type = "Ljung-Box")$p.value
})
plot(lags, pvalues,
     type = "p", pch = 16,
     ylim = c(0, 1),
     xlab = "Lag",
     ylab = "p-value",
     main = "Figure 22. P-values of Ljung-Box Test on Squared Standardised Residuals")
abline(h = 0.05, col = "red", lty = 2)

Figure 22. Ljung-Box p-values for squared GARCH residuals. All points above red line (0.05) confirm no remaining ARCH effects.

# ACF of absolute standardised residuals — further check for remaining structure
acf(abs(residuals(best_garch)),
    na.action = na.omit,
    main = "Figure 23. ACF of Absolute Standardised Residuals - GARCH(1,1)")

Figure 23. ACF of absolute standardised GARCH residuals. No significant autocorrelation confirms model adequacy.

# Normality tests on GARCH standardised residuals
cat("=== Normality Tests on GARCH Residuals ===\n")

## === Normality Tests on GARCH Residuals ===

shapiro.test(na.omit(residuals(best_garch)))

## 
##  Shapiro-Wilk normality test
## 
## data:  na.omit(residuals(best_garch))
## W = 0.99618, p-value = 0.912

jarque.bera.test(na.omit(residuals(best_garch)))

## 
##  Jarque Bera Test
## 
## data:  na.omit(residuals(best_garch))
## X-squared = 0.62537, df = 2, p-value = 0.7315

# Conditional variance plot shows time-varying volatility across study period
plot((fitted(best_garch)[,1])^2, type = "l",
     ylab = "Conditional Variance",
     xlab = "Time",
     main = "Figure 24. Estimated Conditional Variances from GARCH(1,1) - AUD/USD")

Figure 24. Estimated conditional variances from GARCH(1,1). COVID-19 volatility spike clearly visible in early 2020.

---

Section 9 - Forecasting

10-Month ARIMA Forecast

# Generate 10-month ahead point forecast and 80% and 95% prediction intervals
fc_arima <- forecast(best_arima, h = 10, level = c(80, 95))

cat("=== ARIMA(0,1,1) 10-Month Forecast ===\n")

## === ARIMA(0,1,1) 10-Month Forecast ===

print(fc_arima)

##          Point Forecast     Lo 80     Hi 80     Lo 95     Hi 95
## May 2026      0.7081269 0.6755609 0.7406929 0.6583214 0.7579323
## Jun 2026      0.7081269 0.6650956 0.7511581 0.6423162 0.7739375
## Jul 2026      0.7081269 0.6567183 0.7595354 0.6295043 0.7867494
## Aug 2026      0.7081269 0.6495267 0.7667270 0.6185056 0.7977481
## Sep 2026      0.7081269 0.6431259 0.7731278 0.6087165 0.8075373
## Oct 2026      0.7081269 0.6373012 0.7789525 0.5998084 0.8164453
## Nov 2026      0.7081269 0.6319205 0.7843333 0.5915792 0.8246745
## Dec 2026      0.7081269 0.6268953 0.7893584 0.5838940 0.8323598
## Jan 2027      0.7081269 0.6221635 0.7940903 0.5766572 0.8395965
## Feb 2027      0.7081269 0.6176788 0.7985749 0.5697985 0.8464552

# Show last 3 years of history plus 10-month forecast with uncertainty bands
ts_recent <- window(ts_audusd, start = c(2023, 1))

autoplot(ts_recent) +
  autolayer(fc_arima, series = "Forecast", PI = TRUE, alpha = 0.2) +
  geom_line(colour = "#2563EB", linewidth = 1) +
  scale_colour_manual(values = c("Forecast" = "#EF4444")) +
  labs(
    title    = "Figure 25. AUD/USD 10-Month Forecast with 80% and 95% Prediction Intervals",
    subtitle = "ARIMA(0,1,1) + GARCH(1,1)  |  Jan 2023 to Feb 2027",
    x        = "Year",
    y        = "AUD/USD Rate",
    colour   = NULL
  ) +
  theme_minimal(base_size = 13) +
  theme(plot.title      = element_text(face = "bold"),
        legend.position = "bottom")

Figure 25. AUD/USD 10-month forecast (May 2026 to Feb 2027) with 80% and 95% prediction intervals. Widening bands reflect increasing uncertainty at longer horizons.

GARCH Volatility Forecast

# 1-step ahead conditional variance forecast using GARCH(1,1) recursion formula
last_return   <- tail(as.numeric(r_audusd), 1)
last_variance <- tail((fitted(best_garch)[,1])^2, 1)

one_step_var <- garch_coef["a0"] +
                garch_coef["a1"] * last_return^2 +
                garch_coef["b1"] * last_variance

cat("=== GARCH Volatility Forecast ===\n")

## === GARCH Volatility Forecast ===

cat("Last observed conditional variance:", round(last_variance, 4), "\n")

## Last observed conditional variance: 7.7365

cat("1-step ahead forecast of variance :", round(one_step_var, 4), "\n")

## 1-step ahead forecast of variance : 7.9478

cat("1-step ahead forecast of volatility:", round(sqrt(one_step_var), 4), "%\n")

## 1-step ahead forecast of volatility: 2.8192 %

Hold-Out Forecast Accuracy

# Evaluate out-of-sample accuracy by withholding the last 10 months as test set
train_ts <- window(ts_audusd, end   = c(2025, 6))
test_ts  <- window(ts_audusd, start = c(2025, 7))

model_train <- Arima(train_ts, order = arimaorder(best_arima))
fc_test     <- forecast(model_train, h = 10)

cat("=== Forecast Accuracy - Training vs Hold-Out Test Set ===\n")

## === Forecast Accuracy - Training vs Hold-Out Test Set ===

accuracy(fc_test, test_ts)

##                        ME       RMSE        MAE        MPE     MAPE      MASE
## Training set -0.001471658 0.02564035 0.01964137 -0.2461212 2.452361 0.3287290
## Test set      0.021352810 0.03222280 0.02266769  3.0426745 3.245906 0.3793791
##                      ACF1 Theil's U
## Training set 0.0002971361        NA
## Test set     0.5613226313   1.84402

---

Section 10 - Final Combined Forecast Plot

# Final presentation-ready forecast plot combining history and forecast
ts_recent <- window(ts_audusd, start = c(2023, 1))

autoplot(ts_recent) +
  autolayer(fc_arima, series = "Forecast", PI = TRUE, alpha = 0.2) +
  geom_line(colour = "#2563EB", linewidth = 1) +
  scale_colour_manual(values = c("Forecast" = "#EF4444")) +
  labs(
    title    = "Figure 26. AUD/USD - 10-Month Forecast with 95% Confidence Interval",
    subtitle = paste0("ARIMA(", paste(arimaorder(best_arima), collapse = ","),
                      ") + GARCH(1,1)  |  Jan 2023 to Feb 2027"),
    x        = "Year",
    y        = "AUD/USD Rate",
    colour   = NULL
  ) +
  theme_minimal(base_size = 13) +
  theme(plot.title      = element_text(face = "bold"),
        legend.position = "bottom")

Figure 26. Final combined plot showing AUD/USD history from Jan 2023 and 10-month forecast to Feb 2027.

---

Section 11 - Results Summary

# Complete results summary
cat("================================================================\n")

## ================================================================

cat("  RESULTS SUMMARY\n")

##   RESULTS SUMMARY

cat("================================================================\n\n")

## ================================================================

cat("Dataset:\n")

## Dataset:

cat("  Series     : AUD/USD Monthly Exchange Rate\n")

##   Series     : AUD/USD Monthly Exchange Rate

cat("  Source     : Reserve Bank of Australia - F11\n")

##   Source     : Reserve Bank of Australia - F11

cat("  Period     : Jan 2010 to Apr 2026\n")

##   Period     : Jan 2010 to Apr 2026

cat("  n          : 196 observations\n\n")

##   n          : 196 observations

cat("Stationarity:\n")

## Stationarity:

cat("  ADF (original)    p =", round(adf_orig$p.value,  4), "\n")

##   ADF (original)    p = 0.5326

cat("  KPSS (original)   p =", round(kpss_orig$p.value, 4), "\n")

##   KPSS (original)   p = 0.01

cat("  ADF (differenced) p =", round(adf_diff$p.value,  4), "\n")

##   ADF (differenced) p = 0.01

cat("  Decision          : d = 1 applied\n\n")

##   Decision          : d = 1 applied

cat("Seasonality:\n")

## Seasonality:

cat("  No significant seasonal component detected\n")

##   No significant seasonal component detected

cat("  SARIMA excluded; ARIMA + GARCH used instead\n\n")

##   SARIMA excluded; ARIMA + GARCH used instead

cat("Best ARIMA model:\n")

## Best ARIMA model:

cat("  Order     : ARIMA(0,1,1)\n")

##   Order     : ARIMA(0,1,1)

cat("  AIC       :", round(AIC(best_arima), 3), "\n")

##   AIC       : -875.897

cat("  BIC       :", round(BIC(best_arima), 3), "\n")

##   BIC       : -869.351

cat("  Ljung-Box p =", round(lb_test$p.value, 4), "(> 0.05 = PASS)\n\n")

##   Ljung-Box p = 0.3015 (> 0.05 = PASS)

cat("GARCH modelling:\n")

## GARCH modelling:

cat("  McLeod-Li test    : ARCH effects detected (p < 0.05)\n")

##   McLeod-Li test    : ARCH effects detected (p < 0.05)

cat("  Model selected    : GARCH(1,1)\n")

##   Model selected    : GARCH(1,1)

cat("  a1 + b1 =",
    round(garch_coef["a1"] + garch_coef["b1"], 4),
    "(< 1 = mean-reverting)\n\n")

##   a1 + b1 = 0.9591 (< 1 = mean-reverting)

cat("10-Month Point Forecast (ARIMA):\n")

## 10-Month Point Forecast (ARIMA):

print(round(fc_arima$mean, 4))

##         Jan    Feb Mar Apr    May    Jun    Jul    Aug    Sep    Oct    Nov
## 2026                       0.7081 0.7081 0.7081 0.7081 0.7081 0.7081 0.7081
## 2027 0.7081 0.7081                                                         
##         Dec
## 2026 0.7081
## 2027

---

Appendix - AI Tool Usage Declaration

In accordance with the course academic integrity policy, the use of AI tools is declared as follows:

• Tool used: Claude (Anthropic)
• How it was used: Used to assist with R code structure, diagnostic test selection based on course materials (Modules 6, 8, 9), and report writing and formatting.
• Verification of results: All R output, test statistics, model coefficients, and forecasts were verified by running the complete R script independently in RStudio and cross-checked against the course lecture notes.