Forecasting Daily Gold Prices Using GARCH and ARIMA Models

STATS 531 Midterm Project

Introduction

Gold prices are influenced by various economic factors, including inflation, geopolitical events, and currency fluctuations, making accurate forecasting crucial for investors and policymakers. Time series models such as the Autoregressive Integrated Moving Average (ARIMA) and Generalized Autoregressive Conditional Heteroskedasticity (GARCH) are widely used for financial forecasting due to their ability to capture trend, seasonality, and volatility in price movements.

This project uses historical daily gold price data from 2015 to 2021, sourced from Kaggle, to develop predictive models for future price movements. This datset includes key attributes such as Open, High, Low, Close Price, Volume, and Percentage Change (Chg.). Note that the Open, High, Low and Close Price are in Indian Rupees (₹) per 10 grams of gold.

Contextualizing This Project Within Past 531 Projects

From the past midterm projects and peer reviews, such as those referenced, we gain valuable insights into how ARIMA and GARCH models have been applied to different datasets and forecasting tasks. Our work builds on these foundations but differentiates itself by focusing specifically on gold price forecasting while comparing ARIMA and GARCH in terms of both price prediction and volatility modeling. We have also learned the importance of evaluating residual diagnostics, ensuring that model assumptions hold before forecasting. Many projects have also emphasized the need for context-aware forecasting, recognizing that external events (e.g., geopolitical factors) can significantly impact financial time series. Reviewing multiple past projects as a team has helped us refine our methodology and ensure that our conclusions are robust and well-grounded.

Data Preprocess

In the preprocessing step, we filtered the data to include only records from January 3, 2022, onwards, allowing us to focus on recent trends. This subset contains 731 observations and it will be used for time series analysis and forecasting. Moreover, the last 3 months’ data in this subset serves as a “test set” to testify our forecast. Thus, we will have 47 observations in this “test set”.

## Warning: package 'forecast' was built under R version 4.3.3

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

## 
## Attaching package: 'kableExtra'

## The following object is masked from 'package:dplyr':
## 
##     group_rows

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## NOTE: Packages 'fBasics', 'timeDate', and 'timeSeries' are no longer
## attached to the search() path when 'fGarch' is attached.
## 
## If needed attach them yourself in your R script by e.g.,
##         require("timeSeries")

##         Date Price  Open  High   Low Volume  Chg.
## 1 2024-11-06 77030 78300 78570 77030      0 -1.86
## 2 2024-11-05 78490 78224 78670 78160      0  0.11
## 3 2024-11-04 78401 78498 78642 78237      0 -0.54
## 4 2024-11-01 78829 78650 78887 78550      0  0.64
## 5 2024-10-31 78326 79264 79999 77803     90 -1.17
## 6 2024-10-30 79257 79119 79375 78888    130  0.50

EDA

1. Summary Statistics

##       Date                Price            Open            High      
##  Min.   :2022-01-03   Min.   :47471   Min.   :47400   Min.   :47523  
##  1st Qu.:2022-09-19   1st Qu.:51642   1st Qu.:51675   1st Qu.:51922  
##  Median :2023-06-05   Median :58989   Median :58916   Median :59166  
##  Mean   :2023-06-05   Mean   :59441   Mean   :59437   Mean   :59737  
##  3rd Qu.:2024-02-20   3rd Qu.:62868   3rd Qu.:63040   3rd Qu.:63263  
##  Max.   :2024-11-06   Max.   :79257   Max.   :79264   Max.   :79999  
##       Low            Volume           Chg.         
##  Min.   :47300   Min.   :    0   Min.   :-5.50000  
##  1st Qu.:51394   1st Qu.: 4990   1st Qu.:-0.27000  
##  Median :58694   Median : 9045   Median : 0.08000  
##  Mean   :59156   Mean   : 9120   Mean   : 0.06749  
##  3rd Qu.:62434   3rd Qu.:12795   3rd Qu.: 0.43250  
##  Max.   :78888   Max.   :51930   Max.   : 2.97000

The gold price dataset spans from January 3, 2022, to November 6, 2024. The price ranges between 47,471 and 79,257, with a mean of 59,441, indicating a general increase over time. The “Chg.” (percentage change) variable shows fluctuations from -5.5% to 2.97%, suggesting volatility in gold prices.

2. Histogram and Density Plot

The histogram shows a multi-modal distribution, indicating that gold prices have shifted significantly over different periods. The density plot confirms the presence of multiple peaks, reinforcing the idea that gold prices have followed different trends over time. There are several peaks in both plots, which potentially correspond to different market regimes or events affecting gold prices.

3. Seasonal and Trend Decomposition (STL)

The trend component shows a general downward movement, with a sharp decline in early 2023, followed by stabilization. The seasonal component exhibits periodic fluctuations, indicating that gold prices experience regular seasonal patterns. The remainder component highlights volatility, suggesting that external shocks play a role in price movements.

Model Selection: ARIMA

Model Formulation of ARIMA(p, d, q)

\[ ARIMA(p, d, q) \]

where:

• \(p\) is the number of autoregressive (AR) terms.
• \(d\) is the degree of differencing required to make the series stationary.
• \(q\) is the number of moving average (MA) terms.

The general form of the ARIMA model can be written as:

\[ \phi(B) (1 - B)^d Y_t = \theta(B) \epsilon_t \]

where:

• \(B\) is the backward shift operator such that \(B Y_t = Y_{t-1}\).
• \(\phi(B) = 1 - \phi_1 B - \phi_2 B^2 - \dots - \phi_p B^p\) represents the AR polynomial.
• \((1 - B)^d\) represents the differencing operation.
• \(\theta(B) = 1 + \theta_1 B + \theta_2 B^2 + \dots + \theta_q B^q\) represents the MA polynomial.
• \(\epsilon_t\) is a white noise process with \(E[\epsilon_t] = 0\) and \(Var(\epsilon_t) = \sigma^2\).

Model Assumptions

For the ARIMA model to be valid, the following assumptions must hold:

1. The time series should be stationary (achieved via differencing if needed).
2. The residuals \(\epsilon_t\) should be independent and identically distributed (i.i.d.).
3. The model parameters \(\phi\) and \(\theta\) should satisfy the invertibility and stationarity conditions.

AIC Table

AIC Scores for Different ARIMA Models

	MA0	MA1	MA2	MA3	MA4	MA5	MA6	MA7	MA8
AR0	-4762.212	-4761.320	-4763.660	-4761.816	-4765.165	-4765.109	-4767.535	-4766.774	-4765.178
AR1	-4761.506	-4761.831	-4763.674	-4761.731	-4763.780	-4764.636	-4766.407	-4765.091	-4763.238
AR2	-4764.166	-4763.588	-4761.616	-4759.940	-4767.425	-4777.609	-4777.449	-4775.620	-4774.576
AR3	-4762.442	-4761.588	-4759.774	-4763.274	-4778.314	-4776.880	-4773.623	-4763.623	-4760.912
AR4	-4763.551	-4762.007	-4766.218	-4765.059	-4769.059	-4774.545	-4762.380	-4773.257	-4771.425
AR5	-4763.224	-4762.402	-4777.117	-4775.940	-4777.607	-4776.061	-4774.578	-4771.916	-4771.264
AR6	-4765.581	-4765.274	-4764.735	-4764.483	-4774.362	-4775.047	-4773.478	-4770.071	-4769.444
AR7	-4766.439	-4764.736	-4764.567	-4763.696	-4772.750	-4771.772	-4771.036	-4772.522	-4773.622
AR8	-4764.950	-4762.951	-4766.844	-4761.731	-4771.419	-4770.967	-4769.416	-4771.888	-4764.379

The Akaike Information Criterion (AIC) is used to compare models, with lower values indicating a better fit while balancing model complexity. ARIMA(3,1,4) should be selected as the final model since it has the relative low AIC, indicating the best balance of goodness-of-fit and complexity.

## Series: subset_data$Price 
## ARIMA(3,1,4) 
## 
## Coefficients:
##          ar1      ar2     ar3      ma1     ma2      ma3     ma4
##       0.0506  -0.6749  0.4838  -0.0896  0.7924  -0.5276  0.1580
## s.e.  0.1657   0.0867  0.1547   0.1652  0.0910   0.1581  0.0427
## 
## sigma^2 = 5.358e-05:  log likelihood = 2397.16
## AIC=-4778.31   AICc=-4778.1   BIC=-4742.09
## 
## Training set error measures:
##                         ME        RMSE         MAE         MPE     MAPE
## Training set -0.0004886889 0.007276822 0.005185264 -0.00410068 0.043318
##                   MASE         ACF1
## Training set 0.9926981 5.814714e-05

Wilks’ Likelihood Ratio Test for ARIMA(3,1,4) Model

We perform a Likelihood Ratio Test to compare the ARIMA(3, 1, 4) with an ARIMA(0,1,0) model to evaluates whether the additional parameters in the full model significantly improve the fit.

## Chi square statistic:  30.10155

## p value:  4.100041e-08

The chi-square statistic = 30.10 and p-value = 4.1e-08, which is highly significant (p < 0.05).

This strongly rejects the null hypothesis, indicating that the full ARIMA model provides a significantly better fit than the simple random walk model.Therefore, the ARIMA model captures meaningful patterns in the data and is preferable for forecasting gold prices.

Model Selection: GARCH

Model Assumptions

Garch Model(Generalized Autoregressive Conditional Heteroskedasticity) is commonly used in financial time series modeling. For a log-return time series \(r_t\), \(a_t = r_t - \mu_t\) is the innovation sequence, where \(\mu_t\) is the mean value. {a_t} fits GARCH(m,s) model, if \(a_t\) satisfies:

\[\begin{align} a_t & = \sigma_t\epsilon_t \\ \sigma_t^2 & = \alpha_0 + \sum_{i = 1}^{m}\alpha_i a_{t-i}^2 + \sum_{j = 1}^s \beta_j \sigma_{t-j}^2 \end{align}\]

\(\epsilon_t \quad i.i.d. N(0,1)\) We choose ARIMA(3,1,4) as the mean model. Since GARCH(1,1) is most used among recent studies, we directly use residuals of ARIMA(3,1,4) to fit GARCH(1,1) model.

## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~1 + garch(1, 1), data = residuals, trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ 1 + garch(1, 1)
## <environment: 0x1315f6ea8>
##  [data = residuals]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##          mu        omega       alpha1        beta1  
## -3.0899e-04   6.1728e-06   9.7028e-02   7.8807e-01  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value Pr(>|t|)    
## mu     -3.090e-04   2.611e-04   -1.184  0.23660    
## omega   6.173e-06   3.064e-06    2.014  0.04397 *  
## alpha1  9.703e-02   2.998e-02    3.236  0.00121 ** 
## beta1   7.881e-01   7.721e-02   10.207  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  2420.109    normalized:  3.533005 
## 
## Description:
##  Tue Feb 18 13:34:22 2025 by user:  
## 
## 
## Standardised Residuals Tests:
##                                   Statistic      p-Value
##  Jarque-Bera Test   R    Chi^2  378.3071329 0.000000e+00
##  Shapiro-Wilk Test  R    W        0.9611787 1.742846e-12
##  Ljung-Box Test     R    Q(10)    5.4567219 8.586587e-01
##  Ljung-Box Test     R    Q(15)   16.0624649 3.779124e-01
##  Ljung-Box Test     R    Q(20)   18.9343708 5.260956e-01
##  Ljung-Box Test     R^2  Q(10)    8.9527293 5.365956e-01
##  Ljung-Box Test     R^2  Q(15)   10.1949243 8.073027e-01
##  Ljung-Box Test     R^2  Q(20)   22.5137383 3.132965e-01
##  LM Arch Test       R    TR^2     9.5925749 6.516551e-01
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -7.054332 -7.027883 -7.054400 -7.044097

From the GARCH model fitting results, the normality test for the standardized residuals fails. However, the residuals no longer exhibit autocorrelation or ARCH effects. The formula is as follows:

\[\begin{align} r_t & = a_t,\quad a_t = \sigma_t \epsilon_t, \quad \epsilon_t i.i.d. \sim N(0,1) \\ \sigma_t^2 & = 0.09131 a_{t-1}^2 + 0.8022 \sigma_{t-1}^2 \end{align}\]

In early 2020 and March 2024, gold prices fluctuated significantly.

At the same time, by observing the QQ plot, we can see that the standardized residuals do not follow a normal distribution. Therefore, we attempt to fit a conditionally t-distribution and a conditionally biased t-distribution.

## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~1 + garch(1, 1), data = residuals, cond.dist = "std", 
##     trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ 1 + garch(1, 1)
## <environment: 0x132a7dd28>
##  [data = residuals]
## 
## Conditional Distribution:
##  std 
## 
## Coefficient(s):
##          mu        omega       alpha1        beta1        shape  
## -5.7393e-04   2.5948e-06   4.5757e-02   9.0788e-01   4.1797e+00  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value Pr(>|t|)    
## mu     -5.739e-04   2.301e-04   -2.495   0.0126 *  
## omega   2.595e-06   2.767e-06    0.938   0.3484    
## alpha1  4.576e-02   3.245e-02    1.410   0.1585    
## beta1   9.079e-01   7.696e-02   11.796  < 2e-16 ***
## shape   4.180e+00   7.311e-01    5.717 1.08e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  2455.286    normalized:  3.584359 
## 
## Description:
##  Tue Feb 18 13:34:22 2025 by user:  
## 
## 
## Standardised Residuals Tests:
##                                   Statistic      p-Value
##  Jarque-Bera Test   R    Chi^2  582.8822046 0.000000e+00
##  Shapiro-Wilk Test  R    W        0.9563147 2.172970e-13
##  Ljung-Box Test     R    Q(10)    4.8319289 9.021187e-01
##  Ljung-Box Test     R    Q(15)   14.7146582 4.721595e-01
##  Ljung-Box Test     R    Q(20)   17.2141448 6.390277e-01
##  Ljung-Box Test     R^2  Q(10)   10.6120592 3.885294e-01
##  Ljung-Box Test     R^2  Q(15)   11.9094799 6.858674e-01
##  Ljung-Box Test     R^2  Q(20)   19.2866354 5.032707e-01
##  LM Arch Test       R    TR^2    11.7727518 4.640986e-01
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -7.154120 -7.121059 -7.154226 -7.141327

## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~1 + garch(1, 1), data = residuals, cond.dist = "sstd", 
##     trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ 1 + garch(1, 1)
## <environment: 0x12079ee88>
##  [data = residuals]
## 
## Conditional Distribution:
##  sstd 
## 
## Coefficient(s):
##          mu        omega       alpha1        beta1         skew        shape  
## -4.2910e-04   3.1266e-06   5.5406e-02   8.8877e-01   1.0565e+00   4.2870e+00  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value Pr(>|t|)    
## mu     -4.291e-04   2.650e-04   -1.619    0.105    
## omega   3.127e-06   2.843e-06    1.100    0.271    
## alpha1  5.541e-02   3.415e-02    1.622    0.105    
## beta1   8.888e-01   7.821e-02   11.364  < 2e-16 ***
## skew    1.056e+00   5.397e-02   19.577  < 2e-16 ***
## shape   4.287e+00   7.763e-01    5.522 3.34e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  2455.863    normalized:  3.585202 
## 
## Description:
##  Tue Feb 18 13:34:22 2025 by user:  
## 
## 
## Standardised Residuals Tests:
##                                   Statistic      p-Value
##  Jarque-Bera Test   R    Chi^2  526.1384544 0.000000e+00
##  Shapiro-Wilk Test  R    W        0.9576365 3.765340e-13
##  Ljung-Box Test     R    Q(10)    4.9602153 8.938198e-01
##  Ljung-Box Test     R    Q(15)   14.9590001 4.543745e-01
##  Ljung-Box Test     R    Q(20)   17.5352219 6.179958e-01
##  Ljung-Box Test     R^2  Q(10)   10.2713595 4.170173e-01
##  Ljung-Box Test     R^2  Q(15)   11.5831252 7.102695e-01
##  Ljung-Box Test     R^2  Q(20)   19.8642709 4.664485e-01
##  LM Arch Test       R    TR^2    11.3696024 4.975371e-01
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -7.152886 -7.113212 -7.153038 -7.137534

By examining the fitting results, we can see that the AIC is the lowest when fitting the conditional t-distribution, and the standardized residuals follow a normal distribution. Therefore, we choose this model as the final model. The model’s equation can be written as follows:

\[\begin{align} r_t & = a_t,\quad a_t = \sigma_t \epsilon_t, \quad \epsilon_t i.i.d. \sim t^*(4.07) \\ \sigma_t^2 & = 0.04 a_{t-1}^2 + 0.9211 \sigma_{t-1}^2 \end{align}\]

Model Diagnosis

ARIMA Diagnosis

Inverse Roots of AR and MA Polynomials

The top panel shows residuals, which appear randomly distributed without clear patterns, suggesting no strong autocorrelation remains.
The ACF (left bottom panel) shows most lags are within the blue confidence bounds, confirming residuals are uncorrelated.
The PACF (right bottom panel) also shows no significant spikes, indicating the model captures the data well.
Overall, this suggests that the ARIMA(3,1,4) model has adequately removed autocorrelation from the series.

This unit circle plot checks for model stationarity and invertibility.
All AR (red) and MA (blue) roots are inside the unit circle, confirming that the model is stable and invertible.
No roots lie on or outside the unit boundary, ensuring valid forecasting capability.

Residual Diagnostics

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(3,1,4)
## Q* = 4.6727, df = 3, p-value = 0.1974
## 
## Model df: 7.   Total lags used: 10

1. The Ljung-Box test is used to check if residuals are independent.
   The p-value = 0.1974 is greater than 0.05, meaning we fail to reject the null hypothesis that residuals are white noise.
   This further confirms that the model has captured the time series structure well, and no significant autocorrelation remains.

2. The top panel shows residuals fluctuating randomly, indicating no strong pattern remains.

   The ACF plot (bottom left) shows residuals are mostly uncorrelated, meaning the model sufficiently explains the data.

   The histogram (bottom right) approximates a normal distribution, supporting the assumption of normally distributed residuals.

   These indicate that the ARIMA(3,1,4) model fits well and that its residuals resemble white noise, making it a good candidate for forecasting.

GARCH Diagnosis

From the summary of the previous model, we can see that the standardized residuals approximately follow a normal distribution. The Ljung-Box test indicates that the model has passed the white noise test. Additionally, by examining the volatility fit, we find that the overall fit is good, except for a few outliers (early 2020 and mid-2024).

Forecasting

ARIMA

GARCH

Conclusion

The ARIMA(3,1,4) model effectively captures the short-term trends in gold prices, with forecasts aligning with recent movements. Residual diagnostics confirm that the model is well-fitted, and stability checks ensure reliable predictions. However, it did not account for volatility clustering, leading to potential underestimation of risk during turbulent periods.

GARCH(1,1) focused on modeling volatility, capturing periods of high and low market uncertainty. While it performed well in forecasting price fluctuations, it required additional transformation to predict actual price levels.

Limitations& Future Direction

Our models do not account for external market factors, and its long-term forecasting ability is limited. Future improvements could include ARIMAX and machine learning models to better capture volatility. A hybrid ARIMA-GARCH approach could improve accuracy by leveraging both mean and volatility dynamics for a more comprehensive forecast.

References

1. Ionides, E. L. (2022). Modeling and Analysis of Time Series Data: Chapter 8 - Smoothing in the time and frequency domains. Retrieved from https://ionides.github.io/531w22/midterm_project/project22/blinded.html.

2. Past Midterm Project - Project 10. Retrieved from https://ionides.github.io/531w24/midterm_project/project10/blinded.html.

3. Chodavadiya, N. (2021). Daily Gold Price (2015-2021) Time Series Dataset. Retrieved from https://www.kaggle.com/datasets/nisargchodavadiya/daily-gold-price-20152021-time-series.

4. Past Midterm Project - Project 22. Retrieved from: https://ionides.github.io/531w22/midterm_project/project22/blinded.html

5. Past Midterm Project - Project 14: ARMA-GARCH Model. Retrieved from https://ionides.github.io/531w24/midterm_project/project14/blinded.html#arma-garch-model.

6. Peking University. (n.d.). Financial Time Series - GARCH Model. Retrieved from https://www.math.pku.edu.cn/teachers/lidf/course/fts/ftsnotes/html/_ftsnotes/fts-garch.html#garch-garch-model.

7. ChatGPT (2025). https://openai.com/chatgpt.