Time Series Analysis of U.S. Daily Wheat, Corn, and Soybean Futures Prices

I. Introduction

The agricultural commodities market is integral to global trade and economic stability, providing food, feed, and raw materials for a diverse range of industries. However, prices in this sector are highly volatile due to weather-dependent production, geopolitical events, and macroeconomic fluctuations. Futures markets for commodities like wheat, corn, and soybeans provide a vital platform for buyers and sellers to trade standardized contracts for future delivery, primarily facilitating investment opportunities for speculators and hedgers against price volatility.

This document is the full companion code notebook to the written report (report/paper.pdf): it walks through data acquisition and wrangling, descriptive and spectral analysis, and the full modeling pipeline — function-of-time/ARIMA models, univariate GARCH, multivariate (DCC) GARCH, and ARIMAX models — for daily wheat, corn, and soybean futures prices from January 2000 through early December 2024. Each commodity’s futures contract represents 5,000 bushels, with prices quoted in U.S. cents per bushel in the raw data (converted to dollars below).

II. Data Import & Initial Exploration

The code below imports the data for three commodity series with information on prices and volumes of trades since Jan 3, 2000 to Dec 5, 2024. Initial wrangling includes replacing N/A’s values in volumes of each futures with linear approximation.

# Data Imports
wheat <- read.csv("data/wheat_futures.csv")
corn <- read.csv("data/corn_futures.csv")
soy <- read.csv("data/soybean_futures.csv")

wheatvol <- read.csv("data/wheat_volume.csv")
soyvol <- read.csv("data/soybean_volume.csv")

# Fill in N/A's for wheat and soy volume
merged_data <- merge(wheat, wheatvol, by = "Date", all.x = TRUE)
merged_data$Vol.[is.na(merged_data$Vol.)] <- merged_data$Volume[is.na(merged_data$Vol.)]
wheat <- merged_data[, !names(merged_data) %in% "Volume"]

merged_data <- merge(soy, soyvol, by = "Date", all.x = TRUE)
merged_data$Vol.[is.na(merged_data$Vol.)] <- merged_data$Volume[is.na(merged_data$Vol.)]
soy <- merged_data[, !names(merged_data) %in% "Volume"]


# Wrangling
wheat <- wheat %>%
  mutate(Date = as.Date(wheat$Date))

corn <- corn %>%
  mutate(Date = as.Date(corn$Date))

soy <- soy %>%
  mutate(Date = as.Date(soy$Date))


total <- merge(wheat,corn,by="Date")%>%
  rename('Wheat' = 'Price.x',
         'wheat_vol' = 'Vol..x',
         'Corn' = 'Price.y',
         'corn_vol' = 'Vol..y')

total <- merge(total,soy,by="Date")%>%
  rename('Soy' = 'Price',
         'soy_vol' = 'Vol.') %>%
  mutate(wheat_vol = na.approx(wheat_vol))%>%
  mutate(corn_vol = na.approx(corn_vol))%>%
  mutate(soy_vol = na.approx(soy_vol))%>%
  mutate("Month" = month(Date)) %>%
  mutate("Year" = year(Date)) %>%
  mutate(Wheat = Wheat/100)%>%
  mutate(Corn = Corn/100)%>%
  mutate(Soy = Soy/100)

Potential Biases & Assumptions

The data is sourced from the Chicago Board of Trade (CBOT), one of several exchanges (alongside Euronext, KCBT, and MGEX) where these commodities trade, so it represents a sample rather than the full population of futures activity and may not capture dynamics specific to other exchanges. Trading volume figures are also prone to reporting gaps — handled here via linear interpolation — which should be kept in mind when interpreting volume-based results.

The price and volume are plotted below.

# Plot the time series
price_plot <- ggplot(total, aes(x = Date)) +
  geom_line(aes(y = Wheat, color = "Wheat")) +
  geom_line(aes(y = Corn, color = "Corn"))+
  geom_line(aes(y = Soy, color = "Soybeans")) +
  labs(title = "Wheat, Corn, and Soybeans Futures Prices",
       y = "Price ($)",
       x = "Year") +
  scale_color_manual(values = c("Wheat" = "coral", "Corn"= "cornflowerblue","Soybeans"="brown1")) +
  theme_light()

volume_plot <- ggplot(total, aes(x = Date)) +
  geom_line(aes(y = wheat_vol, color = "Wheat"), linewidth = .15) +
  geom_line(aes(y = corn_vol, color = "Corn"), linewidth = .15) +
  geom_line(aes(y = soy_vol, color = "Soybeans"), linewidth = .15) +
  labs(title = "Wheat, Corn, and Soybeans Futures Volumes",
       y = "Volume",
       x = "Year") +
  scale_color_manual(values = c("Wheat" = "coral", "Corn"= "cornflowerblue","Soybeans"="brown1")) +
  scale_y_continuous(labels = scales::label_number(scale_cut = cut_short_scale())) +
  theme_light()

price_plot / volume_plot

As seen below, there does not seem to be any monthly seasonality in the prices of these derivatives.

p <- ggplot(total, aes(x=Month, alpha = Year)) +
  geom_line(aes(y = Wheat, color = "Wheat")) +
  geom_line(aes(y = Corn, color = "Corn"))+
  geom_line(aes(y = Soy, color = "Soybeans"))+
  labs(x = "Month",
       y = "Price ($)",
       title = "Wheat, Corn, and Soybeans Futures Prices by Month") +
  theme_minimal() +
  scale_color_manual(values = c("Wheat" = "coral", "Corn"= "cornflowerblue", "Soybeans" = "brown1"))+
  scale_x_continuous(breaks=c(2,4, 6, 8,10,12)) +
  transition_states(Year)

animate(p, renderer = gifski_renderer())

Now we create time series objects for each commodity’s futures prices and plot them along with corresponding ACF/PACF plots. The most recent 47 trading days (Oct 1 – Dec 5, 2024) are held out as a test set for cross-validation later on.

# Separating Training & Testing data
testdata <- tail(total, 47)
total <- head(total,as.numeric(count(total)-47))

# Constructing xts objects
total.ts <- xts(total, order.by = total$Date)
wheat.ts <- xts(total$Wheat, total$Date)
corn.ts <- xts(total$Corn, total$Date)
soy.ts <- xts(total$Soy, total$Date)


# Original Data
xts_plots(wheat.ts,"Wheat Futures")

xts_plots(corn.ts,"Corn Futures")

xts_plots(soy.ts,"Soybeans Futures")

After plotting the original data, we plot the differenced time series of these derivatives as well.

# Differenced Data
xts_plots(diff(wheat.ts)[-1],"Differenced Wheat Futures")

xts_plots(diff(corn.ts)[-1],"Differenced Corn Futures")

xts_plots(diff(soy.ts)[-1],"Differenced Soy Futures")

When trying to decompose the series, it appears that there is no seasonality that the computer can recognize in either the original or differenced series — decompose()/stl() error out with “time series has no or fewer than 2 periods”, which is expected for a derivatives market rather than a raw seasonal commodity series, so those calls are left commented out below.

# "time series has no or less than 2 periods"

# plot(stl(wheat.ts))
# plot(stl(corn.ts))
# plot(stl(soy.ts))
# plot(stl(wheat.r))
# plot(stl(corn.r))
# plot(stl(soy.r))

Applying Ksmooth, we can try to pick up some overall trend in the data instead.

col.set = c("coral","cornflowerblue","brown1")

plot(na.approx(ts_ts(wheat.ts)), col=col.set[1], lwd=.75, pch=20, ylab="Price ($)", main="Wheat, Corn, and Soybeans Futures - Ksmooth (bandwidth =1)", ylim = c(1.5,18))
lines(ksmooth(time(ts_ts(wheat.ts)), na.approx(ts_ts(wheat.ts)), "normal", bandwidth=1), lwd=2, col=col.set[1])

lines(na.approx(ts_ts(corn.ts)), col=col.set[2], lwd=.75,pch=20)
lines(ksmooth(time(ts_ts(corn.ts)), na.approx(ts_ts(corn.ts)), "normal", bandwidth=1), lwd=2, col=col.set[2])

lines(na.approx(ts_ts(soy.ts)), col=col.set[3], lwd=.75, pch=20)
lines(ksmooth(time(ts_ts(soy.ts)), na.approx(ts_ts(soy.ts)), "normal", bandwidth=1), lwd=2, col=col.set[3])

legend("topleft", col=col.set, lty=1, lwd=2, pch=20, legend=c("Wheat", "Corn","Soybeans"), bg="white", bty = "n")

plot(na.approx(ts_ts(wheat.ts)), col=col.set[1], lwd=.75, pch=20, ylab="Price ($)", main="Wheat, Corn, and Soybeans Futures - Ksmooth (bandwidth =0.5)", ylim = c(1.5,18))
lines(ksmooth(time(ts_ts(wheat.ts)), na.approx(ts_ts(wheat.ts)), "normal", bandwidth=0.5), lwd=2, col=col.set[1])

lines(na.approx(ts_ts(corn.ts)), col=col.set[2], lwd=.75,pch=20)
lines(ksmooth(time(ts_ts(corn.ts)), na.approx(ts_ts(corn.ts)), "normal", bandwidth=0.5), lwd=2, col=col.set[2])

lines(na.approx(ts_ts(soy.ts)), col=col.set[3], lwd=.75, pch=20)
lines(ksmooth(time(ts_ts(soy.ts)), na.approx(ts_ts(soy.ts)), "normal", bandwidth=0.5), lwd=2, col=col.set[3])

legend("topleft", col=col.set, lty=1, lwd=2, pch=20, legend=c("Wheat", "Corn","Soybeans"), bg="white", bty = "n")

III. Spectral Analysis

In this section, we explore the spectrums of each series (differenced for stationarity) and their corresponding daily returns as well.

par(mfrow=c(3,1))
mvspec(diff(wheat.ts)[-1], log = "y", main = "Spectrum of Differenced Wheat Futures (Log)", col = col.set[1], lwd = 2)
mvspec(diff(corn.ts)[-1], log = "y", main = "Spectrum of Differenced Corn Futures (Log)", col = col.set[2], lwd = 2)
mvspec(diff(soy.ts)[-1], log = "y", main = "Spectrum of Differenced Soybeans Futures (Log)", col = col.set[3], lwd = 2)

Bootstrapping the spectral density tells us that there are no significant frequencies.

par(mfrow=c(3,1))
spec1 <- spec.boots(diff(na.approx(ts_ts(wheat.ts))), 1000,kernel("daniell", c(2,2)) )
spec2 <- spec.boots(diff(na.approx(ts_ts(corn.ts))), 1000,kernel("daniell", c(2,2)) )
spec3 <- spec.boots(diff(na.approx(ts_ts(soy.ts))), 1000,kernel("daniell", c(2,2)) )

Daily returns data is produced using the dailyReturn function in quantmod. Each series is stationary, independent, and random. Notably, auto.arima() suggests ARIMA(0,0,0) with zero mean for all of them — i.e., no exploitable structure in the mean of daily returns, which is exactly what motivates modeling their variance with GARCH below.

# Daily returns

# One way to calculate the returns is to use diff(log(data)), which is commented below.
# xts_plots(diff(log(wheat.ts))[-1],"Returns on Wheat Futures",25)

# However, for simplicity, we can also use dailyReturns() from "quantmod".
wheat.r <- dailyReturn(wheat.ts)
xts_plots(wheat.r,"Returns on Wheat Futures",25)

xts_plots(wheat.r^2,"Squared Returns on Wheat Futures",25)

corn.r <- dailyReturn(corn.ts)
xts_plots(corn.r,"Returns on Corn Futures",25)

xts_plots(corn.r^2,"Squared Returns on Corn Futures",25)

soy.r <- dailyReturn(soy.ts)
xts_plots(soy.r,"Returns on Soybeans Futures",25)

xts_plots(soy.r^2,"Squared Returns on Soybeans Futures",25)

# charts.RollingPerformance(R = wheat.r)

summary(auto.arima(wheat.r) )# ARIMA(0,0,0) with zero mean

Series: wheat.r 
ARIMA(0,0,0) with zero mean 

sigma^2 = 0.0004225:  log likelihood = 15452.65
AIC=-30903.3   AICc=-30903.3   BIC=-30896.56

Training set error measures:
                       ME      RMSE        MAE MPE MAPE      MASE        ACF1
Training set 0.0003469623 0.0205549 0.01519133 100  100 0.6987919 0.004276192

summary(auto.arima(corn.r) )# ARIMA(0,0,0) with zero mean

Series: corn.r 
ARIMA(0,0,0) with zero mean 

sigma^2 = 0.0003071:  log likelihood = 16452.73
AIC=-32903.46   AICc=-32903.46   BIC=-32896.71

Training set error measures:
                       ME       RMSE        MAE MPE MAPE      MASE       ACF1
Training set 0.0002679425 0.01752312 0.01248567 100  100 0.6999619 0.01777746

summary(auto.arima(soy.r)) # ARIMA(0,0,0) with zero mean

Series: soy.r 
ARIMA(0,0,0) with zero mean 

sigma^2 = 0.0002496:  log likelihood = 17102.16
AIC=-34202.33   AICc=-34202.33   BIC=-34195.58

Training set error measures:
                       ME       RMSE        MAE MPE MAPE      MASE        ACF1
Training set 0.0002595717 0.01579816 0.01104248 100  100 0.6878481 -0.01606062

par(mfrow=c(3,1))
wheat.r.spec <- mvspec(wheat.r, log = "y", main = "Spectrum of Returns on Wheat Futures (Log)", col = col.set[1], lwd = 1)
corn.r.spec <- mvspec(corn.r, log = "y", main = "Spectrum of Returns on Corn Futures (Log)", col = col.set[2], lwd = 1)
soy.r.spec <- mvspec(soy.r, log = "y", main = "Spectrum of Returns Soybeans Futures (Log)", col = col.set[3], lwd = 1)

IV. ARIMA Models

Wheat Futures

By manually fitting ARIMA models guided by ACF/PACF signals and using auto.arima(), we were not able to find a model with reasonably independent and normal residuals. The “best” model by AIC is ARIMA(1,1,3), suggested by auto.arima().

wheat_auto <- auto.arima(wheat.ts) # ARIMA(1,1,3)
summary(wheat_auto) # AIC=-7404.54 - best

Series: wheat.ts 
ARIMA(1,1,3) 

Coefficients:
         ar1      ma1      ma2      ma3
      0.4723  -0.4298  -0.0275  -0.0485
s.e.  0.1455   0.1456   0.0148   0.0135

sigma^2 = 0.01794:  log likelihood = 3707.27
AIC=-7404.54   AICc=-7404.53   BIC=-7370.83

Training set error measures:
                       ME      RMSE       MAE          MPE     MAPE     MASE
Training set 0.0005738023 0.1338992 0.0858726 -0.007552508 1.520633 1.002313
                     ACF1
Training set 0.0001096424

checkresiduals(wheat_auto)

    Ljung-Box test

data:  Residuals from ARIMA(1,1,3)
Q* = 7.0809, df = 6, p-value = 0.3134

Model df: 4.   Total lags used: 10

acf(wheat_auto$residuals, lag.max = 100)

pacf(wheat_auto$residuals, lag.max = 100)

check_stationarity(wheat_auto$residuals)

Stationarity Test Summary
Test	Test.Statistic	P.Value	Interpretation
ADF Test	-18.991	0.01	Stationary
Phillips-Perron Test	-6274.377	0.01	Stationary
KPSS Test	0.039	0.1	Stationary
ERS Test	-35.844	N/A	Stationary

check_normality_residuals(wheat_auto)

Normality Test Summary
Test	Test.Statistic	P.Value	Interpretation
Anderson-Darling Test	146.469	0	Non-Normal
Jarque-Bera Test	152837.292	0	Non-Normal
Lilliefors Test	0.098	0	Non-Normal
Skewness	0.670	N/A	Right-Skewed
Kurtosis	24.150	N/A	Leptokurtic (Heavy Tails)

wheat_100 <- Arima(wheat.ts, order=c(1,0,0))
summary(wheat_100) #AIC=49867.95

Series: wheat.ts 
ARIMA(1,0,0) with non-zero mean 

Coefficients:
         ar1    mean
      0.9975  5.3236
s.e.  0.0009  0.6314

sigma^2 = 0.01801:  log likelihood = 3692.06
AIC=-7378.13   AICc=-7378.12   BIC=-7357.9

Training set error measures:
                       ME      RMSE        MAE         MPE     MAPE     MASE
Training set 0.0005057067 0.1341914 0.08572483 -0.04281015 1.519925 1.000588
                   ACF1
Training set 0.04524726

checkresiduals(wheat_100)

    Ljung-Box test

data:  Residuals from ARIMA(1,0,0) with non-zero mean
Q* = 42.023, df = 9, p-value = 3.255e-06

Model df: 1.   Total lags used: 10

wheat_010 <- Arima(wheat.ts, order=c(0,1,0))
summary(wheat_010) #AIC=49858.42

Series: wheat.ts 
ARIMA(0,1,0) 

sigma^2 = 0.01803:  log likelihood = 3690.18
AIC=-7378.37   AICc=-7378.37   BIC=-7371.62

Training set error measures:
                       ME     RMSE        MAE          MPE     MAPE     MASE
Training set 0.0005373345 0.134265 0.08566116 -0.007185998 1.515948 0.999845
                   ACF1
Training set 0.04421788

checkresiduals(wheat_010)

    Ljung-Box test

data:  Residuals from ARIMA(0,1,0)
Q* = 42.801, df = 10, p-value = 5.398e-06

Model df: 0.   Total lags used: 10

wheat_110 <- Arima(wheat.ts, order=c(1,1,0))
summary(wheat_110) #AIC=49847.93

Series: wheat.ts 
ARIMA(1,1,0) 

Coefficients:
         ar1
      0.0442
s.e.  0.0126

sigma^2 = 0.018:  log likelihood = 3696.32
AIC=-7388.64   AICc=-7388.63   BIC=-7375.15

Training set error measures:
                       ME      RMSE        MAE          MPE    MAPE     MASE
Training set 0.0005138701 0.1341336 0.08587807 -0.006859185 1.52051 1.002377
                     ACF1
Training set 0.0004391689

checkresiduals(wheat_110)

    Ljung-Box test

data:  Residuals from ARIMA(1,1,0)
Q* = 29.377, df = 9, p-value = 0.0005596

Model df: 1.   Total lags used: 10

wheat_310 <- Arima(wheat.ts, order=c(3,1,0))
summary(wheat_310) #AIC=49834.12

Series: wheat.ts 
ARIMA(3,1,0) 

Coefficients:
         ar1      ar2      ar3
      0.0441  -0.0078  -0.0522
s.e.  0.0126   0.0126   0.0126

sigma^2 = 0.01795:  log likelihood = 3705.2
AIC=-7402.41   AICc=-7402.4   BIC=-7375.44

Training set error measures:
                       ME      RMSE        MAE          MPE     MAPE     MASE
Training set 0.0005465675 0.1339434 0.08590716 -0.007244044 1.521751 1.002716
                     ACF1
Training set -0.001157061

checkresiduals(wheat_310)

    Ljung-Box test

data:  Residuals from ARIMA(3,1,0)
Q* = 11.112, df = 7, p-value = 0.1338

Model df: 3.   Total lags used: 10

wheat.arima <- list("auto.arima() - ARIMA(1,1,3)"= wheat_auto,
                         "ARIMA(1,0,0)" = wheat_100,
                         "ARIMA(0,1,0)" = wheat_010,
                         "ARIMA(1,1,0)" = wheat_110,
                         "ARIMA(3,1,0)" = wheat_310)
arima_summary_table(wheat.arima)

Model Summary
Model	AIC	BIC
auto.arima() - ARIMA(1,1,3)	-7404.54	-7370.83
ARIMA(1,0,0)	-7378.13	-7357.9
ARIMA(0,1,0)	-7378.37	-7371.62
ARIMA(1,1,0)	-7388.64	-7375.15
ARIMA(3,1,0)	-7402.41	-7375.44

Corn Futures

The same thing happens for corn futures — we can’t find a model with desirable residuals, and the best-fit model is again the one suggested by auto.arima(): ARIMA(0,1,3).

corn_auto <- auto.arima(corn.ts) # ARIMA(0,1,3)
summary(corn_auto) # AIC=44553.38 - Best

Series: corn.ts 
ARIMA(0,1,3) 

Coefficients:
         ma1      ma2     ma3
      0.0224  -0.0336  0.0306
s.e.  0.0126   0.0128  0.0124

sigma^2 = 0.006979:  log likelihood = 6665.26
AIC=-13322.51   AICc=-13322.5   BIC=-13295.54

Training set error measures:
                       ME       RMSE        MAE          MPE     MAPE     MASE
Training set 0.0003396775 0.08351438 0.05337484 -0.003958796 1.250147 1.000411
                      ACF1
Training set -0.0003501423

checkresiduals(corn_auto)

    Ljung-Box test

data:  Residuals from ARIMA(0,1,3)
Q* = 6.7363, df = 7, p-value = 0.4568

Model df: 3.   Total lags used: 10

acf(corn_auto$residuals, lag.max = 100)

pacf(corn_auto$residuals, lag.max = 100)

check_stationarity(corn_auto$residuals)

Stationarity Test Summary
Test	Test.Statistic	P.Value	Interpretation
ADF Test	-17.513	0.01	Stationary
Phillips-Perron Test	-6298.253	0.01	Stationary
KPSS Test	0.059	0.1	Stationary
ERS Test	-36.069	N/A	Stationary

check_normality_residuals(corn_auto)

Normality Test Summary
Test	Test.Statistic	P.Value	Interpretation
Anderson-Darling Test	155.914	0	Non-Normal
Jarque-Bera Test	85268.928	0	Non-Normal
Lilliefors Test	0.101	0	Non-Normal
Skewness	-1.140	N/A	Left-Skewed
Kurtosis	17.920	N/A	Leptokurtic (Heavy Tails)

corn_010 <- Arima(corn.ts, order=c(0,1,0))
summary(corn_010) # AIC=44563.78

Series: corn.ts 
ARIMA(0,1,0) 

sigma^2 = 0.006994:  log likelihood = 6657.13
AIC=-13312.25   AICc=-13312.25   BIC=-13305.51

Training set error measures:
                       ME       RMSE        MAE          MPE     MAPE      MASE
Training set 0.0003461912 0.08362282 0.05334472 -0.004067362 1.249022 0.9998467
                   ACF1
Training set 0.01987493

checkresiduals(corn_010)

    Ljung-Box test

data:  Residuals from ARIMA(0,1,0)
Q* = 22.91, df = 10, p-value = 0.01108

Model df: 0.   Total lags used: 10

corn_210 <- Arima(corn.ts, order=c(2,1,0))
summary(corn_210) #AIC=44558.27

Series: corn.ts 
ARIMA(2,1,0) 

Coefficients:
         ar1      ar2
      0.0205  -0.0333
s.e.  0.0126   0.0126

sigma^2 = 0.006986:  log likelihood = 6661.84
AIC=-13317.67   AICc=-13317.67   BIC=-13297.44

Training set error measures:
                       ME       RMSE        MAE          MPE     MAPE     MASE
Training set 0.0003502006 0.08355998 0.05338484 -0.004110672 1.249921 1.000599
                    ACF1
Training set 0.001095628

checkresiduals(corn_210)

    Ljung-Box test

data:  Residuals from ARIMA(2,1,0)
Q* = 13.726, df = 8, p-value = 0.08919

Model df: 2.   Total lags used: 10

corn_213 <- Arima(corn.ts, order=c(2,1,3))
summary(corn_213) #AIC=44554.64

Series: corn.ts 
ARIMA(2,1,3) 

Coefficients:
          ar1      ar2     ma1     ma2     ma3
      -0.4113  -0.4119  0.4337  0.3910  0.0335
s.e.   0.2422   0.2449  0.2419  0.2503  0.0190

sigma^2 = 0.006978:  log likelihood = 6666.58
AIC=-13321.16   AICc=-13321.15   BIC=-13280.71

Training set error measures:
                       ME       RMSE        MAE          MPE     MAPE      MASE
Training set 0.0003398528 0.08349668 0.05335198 -0.003951259 1.249782 0.9999827
                     ACF1
Training set -0.000214407

checkresiduals(corn_213)

    Ljung-Box test

data:  Residuals from ARIMA(2,1,3)
Q* = 4.1716, df = 5, p-value = 0.525

Model df: 5.   Total lags used: 10

corn.arima <- list("auto.arima() - ARIMA(0,1,3)"= corn_auto,
                         "ARIMA(0,1,0)" = corn_010,
                         "ARIMA(2,1,0)" = corn_210,
                         "ARIMA(2,1,3)" = corn_213)

arima_summary_table(corn.arima)

Model Summary
Model	AIC	BIC
auto.arima() - ARIMA(0,1,3)	-13322.51	-13295.54
ARIMA(0,1,0)	-13312.25	-13305.51
ARIMA(2,1,0)	-13317.67	-13297.44
ARIMA(2,1,3)	-13321.16	-13280.71

Soybean Futures

For soybean futures, auto.arima() suggests ARIMA(0,1,0) — effectively a random walk — and we cannot find any better fit.

soy_auto <- auto.arima(soy.ts) # ARIMA(0,1,0)
summary(soy_auto) # AIC=44553.38 - Best

Series: soy.ts 
ARIMA(0,1,0) 

sigma^2 = 0.03026:  log likelihood = 2067.51
AIC=-4133.02   AICc=-4133.02   BIC=-4126.28

Training set error measures:
                       ME      RMSE       MAE          MPE     MAPE     MASE
Training set 0.0009589221 0.1739521 0.1116886 0.0007420328 1.106765 0.999847
                    ACF1
Training set -0.01909734

checkresiduals(soy_auto)

    Ljung-Box test

data:  Residuals from ARIMA(0,1,0)
Q* = 25.398, df = 10, p-value = 0.00464

Model df: 0.   Total lags used: 10

acf(soy_auto$residuals, lag.max = 100)

pacf(soy_auto$residuals, lag.max = 100)

check_stationarity(soy_auto$residuals)

Stationarity Test Summary
Test	Test.Statistic	P.Value	Interpretation
ADF Test	-18.627	0.01	Stationary
Phillips-Perron Test	-6302.494	0.01	Stationary
KPSS Test	0.068	0.1	Stationary
ERS Test	-37.369	N/A	Stationary

check_normality_residuals(soy_auto)

Normality Test Summary
Test	Test.Statistic	P.Value	Interpretation
Anderson-Darling Test	143.728	0	Non-Normal
Jarque-Bera Test	275226.757	0	Non-Normal
Lilliefors Test	0.104	0	Non-Normal
Skewness	-1.080	N/A	Left-Skewed
Kurtosis	32.380	N/A	Leptokurtic (Heavy Tails)

soy_110 <- Arima(soy.ts, order=c(0,1,1))
summary(soy_110) #AIC=44554.64

Series: soy.ts 
ARIMA(0,1,1) 

Coefficients:
          ma1
      -0.0193
s.e.   0.0127

sigma^2 = 0.03026:  log likelihood = 2068.67
AIC=-4133.33   AICc=-4133.33   BIC=-4119.85

Training set error measures:
                       ME    RMSE       MAE          MPE     MAPE      MASE
Training set 0.0009780841 0.17392 0.1116465 0.0007637643 1.106114 0.9994705
                     ACF1
Training set 0.0001024806

checkresiduals(soy_110)

    Ljung-Box test

data:  Residuals from ARIMA(0,1,1)
Q* = 23.742, df = 9, p-value = 0.004729

Model df: 1.   Total lags used: 10

V. GARCH Models

Wheat Futures

We fit ARMA(0,0)+GARCH(1,1) and ARMA(1,1)+GARCH(1,1). Running garch_compare() to search over conditional error distributions points to a skewed Student-t (“sstd”). This differs from the ARMA(1,3)+GARCH(1,1) suggested by garchAuto()’s own AIC-based search, but the simpler ARMA(0,0)+GARCH(1,1) specification is preferred here since it aligns with the returns analysis above (auto.arima() suggesting ARIMA(0,0,0) for the mean).

# Returns of Wheat

wheat.garch_auto = garchAuto(wheat.r, cores=8, trace=TRUE)
summary(wheat.garch_auto)

   35.71326 0.001944128
 Ljung-Box Test     R^2  Q(20)   42.52310 0.002362458
 LM Arch Test       R    TR^2    27.50155 0.006539122

Information Criterion Statistics:
      AIC       BIC       SIC      HQIC 
-5.097382 -5.085546 -5.097388 -5.093280 

mod_specify = ugarchspec(mean.model = list(armaOrder=c(0,0)),variance.model = list(model = "sGARCH",garchOrder=c(1,1)),distribution.model='sstd')

wheat.garch1 = ugarchfit(data = wheat.r, spec = mod_specify, out.sample = 20)
wheat.garch1

*---------------------------------*
*          GARCH Model Fit        *
*---------------------------------*

Conditional Variance Dynamics   
-----------------------------------
GARCH Model : sGARCH(1,1)
Mean Model  : ARFIMA(0,0,0)
Distribution    : sstd 

Optimal Parameters
------------------------------------
        Estimate  Std. Error  t value Pr(>|t|)
mu      0.000267    0.000228   1.1724 0.241017
omega   0.000005    0.000003   1.7507 0.079999
alpha1  0.044278    0.008461   5.2332 0.000000
beta1   0.944618    0.011798  80.0690 0.000000
skew    1.170023    0.022124  52.8855 0.000000
shape   8.280267    0.783354  10.5703 0.000000

Robust Standard Errors:
        Estimate  Std. Error  t value Pr(>|t|)
mu      0.000267    0.000221   1.2065  0.22763
omega   0.000005    0.000010   0.4578  0.64709
alpha1  0.044278    0.027800   1.5927  0.11122
beta1   0.944618    0.042515  22.2184  0.00000
skew    1.170023    0.026454  44.2278  0.00000
shape   8.280267    0.814105  10.1710  0.00000

LogLikelihood : 15939.36 

Information Criteria
------------------------------------
                    
Akaike       -5.1011
Bayes        -5.0947
Shibata      -5.1011
Hannan-Quinn -5.0989

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                        statistic p-value
Lag[1]                     0.6311  0.4269
Lag[2*(p+q)+(p+q)-1][2]    0.7094  0.6029
Lag[4*(p+q)+(p+q)-1][5]    0.9320  0.8749
d.o.f=0
H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic  p-value
Lag[1]                      1.485 0.222972
Lag[2*(p+q)+(p+q)-1][5]    10.260 0.007892
Lag[4*(p+q)+(p+q)-1][9]    15.916 0.002140
d.o.f=2

Weighted ARCH LM Tests
------------------------------------
            Statistic Shape Scale  P-Value
ARCH Lag[3]     5.303 0.500 2.000 0.021290
ARCH Lag[5]    11.672 1.440 1.667 0.002558
ARCH Lag[7]    13.199 2.315 1.543 0.003188

Nyblom stability test
------------------------------------
Joint Statistic:  16.9177
Individual Statistics:             
mu     0.1007
omega  1.9032
alpha1 0.1465
beta1  0.1219
skew   0.1857
shape  0.2397

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:         1.49 1.68 2.12
Individual Statistic:    0.35 0.47 0.75

Sign Bias Test
------------------------------------
                   t-value    prob sig
Sign Bias           0.1320 0.89496    
Negative Sign Bias  0.2979 0.76579    
Positive Sign Bias  2.2042 0.02754  **
Joint Effect        9.4444 0.02393  **


Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
  group statistic p-value(g-1)
1    20     36.67     0.008718
2    30     46.82     0.019419
3    40     70.28     0.001566
4    50     66.25     0.050792


Elapsed time : 0.258523 

checkresiduals(wheat.garch1@fit)

    Ljung-Box test

data:  Residuals
Q* = 3.7335, df = 10, p-value = 0.9586

Model df: 0.   Total lags used: 10

acf2(wheat.garch1@fit$residuals^2)

     [,1] [,2] [,3] [,4] [,5]  [,6] [,7]  [,8] [,9] [,10] [,11] [,12] [,13]
ACF  0.15 0.15 0.22 0.17 0.18  0.08 0.10  0.08 0.12  0.08  0.05  0.06  0.09
PACF 0.15 0.13 0.19 0.10 0.11 -0.01 0.02 -0.01 0.07  0.01  0.00  0.00  0.04
     [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
ACF   0.07  0.05  0.02  0.05  0.04  0.03  0.03  0.02  0.05  0.03  0.03  0.05
PACF  0.01  0.01 -0.03  0.01  0.00  0.01  0.00  0.00  0.02  0.01  0.01  0.03
     [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37]
ACF   0.03  0.05  0.03  0.03  0.04  0.02  0.03  0.01  0.03  0.02  0.02  0.03
PACF -0.01  0.02 -0.01  0.01  0.01  0.00  0.00 -0.01  0.00  0.01  0.00  0.02
     [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48] [,49]
ACF   0.02  0.02  0.03  0.03  0.01  0.02  0.04  0.04  0.03  0.05  0.01  0.03
PACF  0.01  0.00  0.00  0.02 -0.01  0.01  0.02  0.02  0.01  0.03 -0.02  0.00
     [,50] [,51] [,52] [,53] [,54] [,55] [,56] [,57] [,58] [,59] [,60] [,61]
ACF   0.05  0.03  0.05  0.04  0.04  0.03  0.04  0.03  0.03  0.03  0.03  0.01
PACF  0.02  0.01  0.02  0.01  0.01 -0.01  0.00  0.00  0.01  0.00  0.00 -0.02
     [,62] [,63] [,64] [,65] [,66] [,67] [,68] [,69] [,70] [,71] [,72] [,73]
ACF   0.02  0.06  0.03  0.02  0.04  0.05  0.03  0.01  0.02  0.02  0.04  0.03
PACF  0.00  0.04  0.01  0.00  0.01  0.02 -0.01 -0.02  0.00  0.00  0.02  0.01
     [,74] [,75] [,76] [,77] [,78] [,79] [,80] [,81] [,82] [,83] [,84] [,85]
ACF   0.05  0.02  0.05  0.04  0.02  0.03  0.02  0.03  0.02  0.03  0.04  0.06
PACF  0.03 -0.01  0.02  0.00 -0.01  0.00 -0.01  0.01  0.00  0.01  0.03  0.04
     [,86] [,87] [,88] [,89] [,90]
ACF   0.04  0.03  0.04  0.04  0.02
PACF  0.01  0.00  0.00  0.00 -0.01

plot(wheat.garch1, which=1)

mod_specify = ugarchspec(mean.model = list(armaOrder=c(1,1)),variance.model = list(model = "sGARCH",garchOrder=c(1,1)),distribution.model='sstd')

wheat.garch2 = ugarchfit(data = wheat.r, spec = mod_specify, out.sample = 20)
wheat.garch2

*---------------------------------*
*          GARCH Model Fit        *
*---------------------------------*

Conditional Variance Dynamics   
-----------------------------------
GARCH Model : sGARCH(1,1)
Mean Model  : ARFIMA(1,0,1)
Distribution    : sstd 

Convergence Problem:
Solver Message: 

garch_compare(wheat.r, arma_order = c(0, 0), garch_order = c(1, 1))

Warning - Persistence: 1.152427 

GARCH Model Comparison
Model	AIC	BIC
ARMA(0,0)-GARCH(1,1) with std	-5.092	-5.086
ARMA(0,0)-GARCH(1,1) with snorm	-5.075	-5.07
ARMA(0,0)-GARCH(1,1) with ged	-5.084	-5.078
ARMA(0,0)-GARCH(1,1) with sged	-5.096	-5.09
ARMA(0,0)-GARCH(1,1) with sstd	-5.103	-5.096
ARMA(0,0)-GARCH(1,1) with snig	-5.102	-5.096
ARMA(0,0)-GARCH(1,1) with QMLE	-5.06	-5.056

Residuals are largely independent and normal.

checkresiduals(wheat.garch1@fit)

    Ljung-Box test

data:  Residuals
Q* = 3.7335, df = 10, p-value = 0.9586

Model df: 0.   Total lags used: 10

acf(wheat.garch1@fit$residuals, lag.max=100)

pacf(wheat.garch1@fit$residuals, lag.max=100)

This plot shows the relationship between the squared residuals and the model’s fitted conditional variance — the two tracking each other is a sanity check that the GARCH variance equation is picking up the volatility clustering visible in the data.

plot((wheat.garch1@fit$residuals)^2, type = "l",ylab="Squared Residuals")
lines(wheat.garch1@fit$var , col = "green")

Corn and Soybean Futures

We repeat the same process for returns on corn and soybean futures. The best model for each series is again ARMA(0,0)+GARCH(1,1).

# Returns of Corn

corn_auto = garchAuto(corn.r, cores=8, trace=TRUE)
summary(corn_auto)

 R^2  Q(10)     7.240327 0.7025785
 Ljung-Box Test     R^2  Q(15)    10.758165 0.7695417
 Ljung-Box Test     R^2  Q(20)    18.182881 0.5753623
 LM Arch Test       R    TR^2     10.036841 0.6127285

Information Criterion Statistics:
      AIC       BIC       SIC      HQIC 
-5.473799 -5.458736 -5.473809 -5.468579 

garch_compare(corn.r, arma_order = c(1, 1), garch_order = c(1, 1))

GARCH Model Comparison
Model	AIC	BIC
ARMA(1,1)-GARCH(1,1) with std	-5.481	-5.473
ARMA(1,1)-GARCH(1,1) with snorm	-5.406	-5.398
ARMA(1,1)-GARCH(1,1) with ged	-5.47	-5.462
ARMA(1,1)-GARCH(1,1) with sged	-5.471	-5.462
ARMA(1,1)-GARCH(1,1) with sstd	-5.482	-5.473
ARMA(1,1)-GARCH(1,1) with QMLE	-5.405	-5.398

garch_compare(corn.r, arma_order = c(0, 0), garch_order = c(1, 1))

GARCH Model Comparison
Model	AIC	BIC
ARMA(0,0)-GARCH(1,1) with std	-5.481	-5.476
ARMA(0,0)-GARCH(1,1) with snorm	-5.406	-5.401
ARMA(0,0)-GARCH(1,1) with ged	-5.47	-5.465
ARMA(0,0)-GARCH(1,1) with sged	-5.471	-5.464
ARMA(0,0)-GARCH(1,1) with sstd	-5.482	-5.475
ARMA(0,0)-GARCH(1,1) with QMLE	-5.405	-5.4

mod_specify = ugarchspec(mean.model = list(armaOrder=c(0,0)),variance.model = list(model = "sGARCH",garchOrder=c(1,1)),distribution.model='std')

corn.garch1 = ugarchfit(data = corn.r, spec = mod_specify, out.sample = 20)
corn.garch1

*---------------------------------*
*          GARCH Model Fit        *
*---------------------------------*

Conditional Variance Dynamics   
-----------------------------------
GARCH Model : sGARCH(1,1)
Mean Model  : ARFIMA(0,0,0)
Distribution    : std 

Optimal Parameters
------------------------------------
        Estimate  Std. Error  t value Pr(>|t|)
mu      0.000037    0.000171  0.21433 0.830289
omega   0.000004    0.000002  2.09812 0.035894
alpha1  0.067258    0.010450  6.43607 0.000000
beta1   0.920684    0.013055 70.52169 0.000000
shape   5.340883    0.352992 15.13034 0.000000

Robust Standard Errors:
        Estimate  Std. Error  t value Pr(>|t|)
mu      0.000037    0.000165  0.22231 0.824075
omega   0.000004    0.000007  0.65450 0.512788
alpha1  0.067258    0.025224  2.66641 0.007667
beta1   0.920684    0.036939 24.92423 0.000000
shape   5.340883    0.397714 13.42894 0.000000

LogLikelihood : 17130.95 

Information Criteria
------------------------------------
                    
Akaike       -5.4829
Bayes        -5.4775
Shibata      -5.4829
Hannan-Quinn -5.4811

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                        statistic p-value
Lag[1]                      1.750  0.1859
Lag[2*(p+q)+(p+q)-1][2]     2.771  0.1619
Lag[4*(p+q)+(p+q)-1][5]     3.927  0.2633
d.o.f=0
H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                     0.4937  0.4823
Lag[2*(p+q)+(p+q)-1][5]    1.4426  0.7540
Lag[4*(p+q)+(p+q)-1][9]    3.4916  0.6758
d.o.f=2

Weighted ARCH LM Tests
------------------------------------
            Statistic Shape Scale P-Value
ARCH Lag[3]   0.07015 0.500 2.000  0.7911
ARCH Lag[5]   1.22112 1.440 1.667  0.6686
ARCH Lag[7]   2.51592 2.315 1.543  0.6092

Nyblom stability test
------------------------------------
Joint Statistic:  13.2548
Individual Statistics:             
mu     0.1301
omega  1.3874
alpha1 0.4360
beta1  0.4841
shape  0.4789

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:         1.28 1.47 1.88
Individual Statistic:    0.35 0.47 0.75

Sign Bias Test
------------------------------------
                   t-value    prob sig
Sign Bias           1.3664 0.17186    
Negative Sign Bias  0.6427 0.52045    
Positive Sign Bias  0.6773 0.49824    
Joint Effect        9.2907 0.02567  **


Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
  group statistic p-value(g-1)
1    20     43.02    1.287e-03
2    30     70.19    2.858e-05
3    40     96.84    7.993e-07
4    50    144.16    2.644e-11


Elapsed time : 0.138484 

mod_specify = ugarchspec(mean.model = list(armaOrder=c(1,1)),variance.model = list(model = "sGARCH",garchOrder=c(1,1)),distribution.model='std')

corn.garch2 = ugarchfit(data = corn.r, spec = mod_specify, out.sample = 20)
corn.garch2

*---------------------------------*
*          GARCH Model Fit        *
*---------------------------------*

Conditional Variance Dynamics   
-----------------------------------
GARCH Model : sGARCH(1,1)
Mean Model  : ARFIMA(1,0,1)
Distribution    : std 

Optimal Parameters
------------------------------------
        Estimate  Std. Error    t value Pr(>|t|)
mu      0.000030    0.000172    0.17272 0.862869
ar1    -0.933477    0.002313 -403.49675 0.000000
ma1     0.941863    0.004579  205.69590 0.000000
omega   0.000004    0.000002    2.24423 0.024817
alpha1  0.067259    0.009612    6.99716 0.000000
beta1   0.920725    0.011856   77.65630 0.000000
shape   5.339638    0.352715   15.13867 0.000000

Robust Standard Errors:
        Estimate  Std. Error    t value Pr(>|t|)
mu      0.000030    0.000166    0.17937 0.857644
ar1    -0.933477    0.007988 -116.86349 0.000000
ma1     0.941863    0.005896  159.75731 0.000000
omega   0.000004    0.000006    0.74365 0.457088
alpha1  0.067259    0.020643    3.25813 0.001121
beta1   0.920725    0.030735   29.95677 0.000000
shape   5.339638    0.394279   13.54278 0.000000

LogLikelihood : 17132.77 

Information Criteria
------------------------------------
                    
Akaike       -5.4829
Bayes        -5.4753
Shibata      -5.4829
Hannan-Quinn -5.4803

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                        statistic p-value
Lag[1]                     0.6033  0.4373
Lag[2*(p+q)+(p+q)-1][5]    1.4812  0.9983
Lag[4*(p+q)+(p+q)-1][9]    3.7728  0.7430
d.o.f=2
H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                      0.459  0.4981
Lag[2*(p+q)+(p+q)-1][5]     1.405  0.7633
Lag[4*(p+q)+(p+q)-1][9]     3.441  0.6845
d.o.f=2

Weighted ARCH LM Tests
------------------------------------
            Statistic Shape Scale P-Value
ARCH Lag[3]     0.098 0.500 2.000  0.7542
ARCH Lag[5]     1.217 1.440 1.667  0.6698
ARCH Lag[7]     2.491 2.315 1.543  0.6142

Nyblom stability test
------------------------------------
Joint Statistic:  13.9429
Individual Statistics:             
mu     0.1262
ar1    0.0993
ma1    0.1110
omega  1.4219
alpha1 0.4316
beta1  0.4851
shape  0.4675

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:         1.69 1.9 2.35
Individual Statistic:    0.35 0.47 0.75

Sign Bias Test
------------------------------------
                   t-value    prob sig
Sign Bias           1.3365 0.18144    
Negative Sign Bias  0.6268 0.53083    
Positive Sign Bias  0.7613 0.44652    
Joint Effect        9.4040 0.02438  **


Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
  group statistic p-value(g-1)
1    20     18.40       0.4960
2    30     27.20       0.5611
3    40     29.32       0.8697
4    50     60.53       0.1250


Elapsed time : 0.2587819 

plot(corn.garch2, which=1)

mod_specify = ugarchspec(mean.model = list(armaOrder=c(1,1)),variance.model = list(model = "eGARCH",garchOrder=c(0,1)),distribution.model='std')

corn.egarch = ugarchfit(data = corn.r, spec = mod_specify, out.sample = 20)
corn.egarch

*---------------------------------*
*          GARCH Model Fit        *
*---------------------------------*

Conditional Variance Dynamics   
-----------------------------------
GARCH Model : eGARCH(0,1)
Mean Model  : ARFIMA(1,0,1)
Distribution    : std 

Optimal Parameters
------------------------------------
       Estimate  Std. Error     t value Pr(>|t|)
mu     0.000001    0.000187  4.9350e-03  0.99606
ar1   -0.939532    0.016638 -5.6468e+01  0.00000
ma1    0.947168    0.014740  6.4261e+01  0.00000
omega -0.010306    0.000034 -2.9876e+02  0.00000
beta1  0.998703    0.000001  8.7862e+05  0.00000
shape  3.566747    0.010353  3.4453e+02  0.00000

Robust Standard Errors:
       Estimate  Std. Error     t value Pr(>|t|)
mu     0.000001    0.000180  5.1140e-03  0.99592
ar1   -0.939532    0.005535 -1.6975e+02  0.00000
ma1    0.947168    0.007234  1.3094e+02  0.00000
omega -0.010306    0.000096 -1.0728e+02  0.00000
beta1  0.998703    0.000002  6.0320e+05  0.00000
shape  3.566747    0.089459  3.9870e+01  0.00000

LogLikelihood : 16863.38 

Information Criteria
------------------------------------
                    
Akaike       -5.3970
Bayes        -5.3905
Shibata      -5.3970
Hannan-Quinn -5.3947

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                        statistic p-value
Lag[1]                     0.5739  0.4487
Lag[2*(p+q)+(p+q)-1][5]    3.1610  0.3744
Lag[4*(p+q)+(p+q)-1][9]    5.4089  0.3653
d.o.f=2
H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic   p-value
Lag[1]                      45.34 1.659e-11
Lag[2*(p+q)+(p+q)-1][2]     70.13 0.000e+00
Lag[4*(p+q)+(p+q)-1][5]    127.29 0.000e+00
d.o.f=1

Weighted ARCH LM Tests
------------------------------------
            Statistic Shape Scale   P-Value
ARCH Lag[2]     49.55 0.500 2.000 1.932e-12
ARCH Lag[4]     89.82 1.397 1.611 0.000e+00
ARCH Lag[6]    123.74 2.222 1.500 0.000e+00

Nyblom stability test
------------------------------------
Joint Statistic:  16.6769
Individual Statistics:            
mu    0.1068
ar1   0.1136
ma1   0.1097
omega 5.7932
beta1 5.7766
shape 3.7512

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:         1.49 1.68 2.12
Individual Statistic:    0.35 0.47 0.75

Sign Bias Test
------------------------------------
                   t-value      prob sig
Sign Bias           0.4078 6.834e-01    
Negative Sign Bias  7.1859 7.455e-13 ***
Positive Sign Bias  3.8513 1.187e-04 ***
Joint Effect       69.5696 5.277e-15 ***


Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
  group statistic p-value(g-1)
1    20     31.63      0.03436
2    30     34.46      0.22288
3    40     42.98      0.30478
4    50     53.58      0.30283


Elapsed time : 0.4054549 

checkresiduals(corn.garch1@fit)

    Ljung-Box test

data:  Residuals
Q* = 10.485, df = 10, p-value = 0.399

Model df: 0.   Total lags used: 10

acf(corn.garch1@fit$residuals, lag.max=100)

pacf(corn.garch1@fit$residuals, lag.max=100)

soy.garch_auto = garchAuto(soy.r, cores=8, trace=TRUE)
summary(soy.garch_auto)


mod_specify = ugarchspec(mean.model = list(armaOrder=c(0,0)),variance.model = list(model = "sGARCH",garchOrder=c(1,1)),distribution.model='std')

soy.garch1 = ugarchfit(data = soy.r, spec = mod_specify, out.sample = 30)
soy.garch1

     19.84 2.202e-04
d.o.f=2

Weighted ARCH LM Tests
------------------------------------
            Statistic Shape Scale P-Value
ARCH Lag[3] 0.0003901 0.500 2.000  0.9842
ARCH Lag[5] 2.4388767 1.440 1.667  0.3822
ARCH Lag[7] 2.9599188 2.315 1.543  0.5212

Nyblom stability test
------------------------------------
Joint Statistic:  79.0575
Individual Statistics:             
mu     0.3279
omega  8.5845
alpha1 0.5573
beta1  0.7205
shape  0.5777

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:         1.28 1.47 1.88
Individual Statistic:    0.35 0.47 0.75

Sign Bias Test
------------------------------------
                   t-value      prob sig
Sign Bias            1.843 0.0653079   *
Negative Sign Bias   2.089 0.0367240  **
Positive Sign Bias   3.069 0.0021551 ***
Joint Effect        16.982 0.0007129 ***


Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
  group statistic p-value(g-1)
1    20     21.65       0.3019
2    30     32.89       0.2820
3    40     35.36       0.6369
4    50     48.47       0.4947


Elapsed time : 0.1509089 

checkresiduals(soy.garch1@fit)

    Ljung-Box test

data:  Residuals
Q* = 20.965, df = 10, p-value = 0.02134

Model df: 0.   Total lags used: 10

checkresiduals(soy.garch1@fit)

    Ljung-Box test

data:  Residuals
Q* = 20.965, df = 10, p-value = 0.02134

Model df: 0.   Total lags used: 10

acf(soy.garch1@fit$residuals, lag.max=100)

pacf(soy.garch1@fit$residuals, lag.max=100)

VI. GARCH Forecasting & Cross-Validation

Using the GARCH models developed above, we forecast into the test window (Oct 1 – Dec 5, 2024) held out earlier. That window is 47 trading days, but ts_ts() below regularizes each irregular trading-day return series onto a full calendar-day scale (inserting NAs for weekends/holidays) — so the actual comparison series ends up 66 points long, not 47. The forecast horizon (test_horizon) is therefore set to the calendar-day span of the test window, not its row count, so the predicted and actual series line up one-to-one for the R² calculation below.

test_horizon <- as.numeric(diff(range(testdata$Date))) + 1

# Forecast Wheat
wheat.for <- ugarchforecast(wheat.garch1, n.ahead = test_horizon)
wheat.test <- ts_ts(dailyReturn(xts(testdata$Wheat,testdata$Date)))

wheat.for.returns <- wheat.for@forecast$seriesFor
wheat_r_t <- c(tail(ts_ts(wheat.r),150),rep(NA,test_horizon))
wheat_r_f <- c(rep(NA,150),wheat.for.returns)
wheat_r_a <- c(rep(NA,150),wheat.test)

# Forecast Corn
corn.for <- ugarchforecast(corn.garch1, n.ahead = test_horizon)
corn.test <-  ts_ts(dailyReturn(xts(testdata$Corn,testdata$Date)))

corn.for.returns <- corn.for@forecast$seriesFor
corn_r_t <- c(tail(ts_ts(corn.r),150),rep(NA,test_horizon))
corn_r_f <- c(rep(NA,150),corn.for.returns)
corn_r_a <- c(rep(NA,150),corn.test)


# Forecast Soybeans
soy.for <- ugarchforecast(soy.garch1, n.ahead = test_horizon)
soy.test <-  ts_ts(dailyReturn(xts(testdata$Soy,testdata$Date)))

soy.for.returns <- soy.for@forecast$seriesFor
soy_r_t <- c(tail(ts_ts(soy.r),150),rep(NA,test_horizon))
soy_r_f <- c(rep(NA,150),soy.for.returns)
soy_r_a <- c(rep(NA,150),soy.test)



# Plotting
par(mfrow=c(3,1))
plot(wheat_r_t, type = "l", main = "Forecasted Wheat Futures Returns")
lines(wheat_r_f, col = "orange")
lines(wheat_r_a, col = "green")

plot(corn_r_t, type = "l", main = "Forecasted Corn Futures Returns")
lines(corn_r_f, col = "orange")
lines(corn_r_a, col = "green")

plot(soy_r_t, type = "l", main = "Forecasted Soybeans Futures Returns")
lines(soy_r_f, col = "orange")
lines(soy_r_a, col = "green")

Unfortunately, all three models yield negative R-Squared on this held-out window, meaning the forecasted values perform worse than simply predicting the mean return.

# Calculate R-squared
wheat_forecast_df <- na.omit(data.frame(
  predicted_values = wheat_r_f,
  actual_values = wheat_r_a))

calculate_r2(
  predicted_values = wheat_forecast_df$predicted_values,
  actual_values = wheat_forecast_df$actual_values)

R-squared of the model is: -0.0123 

corn_forecast_df <- na.omit(data.frame(
  predicted_values = corn_r_f,
  actual_values = corn_r_a))

calculate_r2(
  predicted_values = corn_forecast_df$predicted_values,
  actual_values = corn_forecast_df$actual_values)

R-squared of the model is: -9e-04 

soy_forecast_df <- na.omit(data.frame(
  predicted_values = soy_r_f,
  actual_values = soy_r_a))

calculate_r2(
  predicted_values = soy_forecast_df$predicted_values,
  actual_values = soy_forecast_df$actual_values)

R-squared of the model is: -0.0329 

VII. Multivariate GARCH

A DCC (dynamic conditional correlation) multivariate GARCH model is fit across all three return series to estimate how the correlation between each pair of commodities evolves over time, rather than assuming it is constant.

rX <- data.frame(wheat.r, corn.r, soy.r)
names(rX)[1] <- "wheat.r"
names(rX)[2] <- "corn.r"
names(rX)[3] <- "soy.r"

# DCC (MVN)
uspec.n = multispec(replicate(3, ugarchspec(mean.model = list(armaOrder = c(0,0)))))
multf = multifit(uspec.n, rX)
spec1 = dccspec(uspec = uspec.n, dccOrder = c(1, 1), distribution = 'mvnorm')

fit1 = dccfit(spec1, data = rX, fit.control = list(eval.se = TRUE), fit = multf)
cov1 = rcov(fit1)  # extracts the covariance matrix
cor1 = rcor(fit1)  # extracts the correlation matrix

par(mfrow=c(3,1))
plot(as.xts(cor1[1,2,]),main="Wheat and Corn")
plot(as.xts(cor1[1,3,]),main="Wheat and Soybeans")
plot(as.xts(cor1[2,3,]),main="Corn and Soybeans")

The forecast correlation is shown below, using the same test_horizon as the GARCH forecasts above for consistency.

dccf1 <- dccforecast(fit1, n.ahead = test_horizon)
Rf <- dccf1@mforecast$R
corf_WC <- Rf[[1]][1,2,]
corf_WS <- Rf[[1]][1,3,]
corf_CS <- Rf[[1]][2,3,]




par(mfrow=c(3,1))
c_WC <- c(tail(cor1[1,2,],150),rep(NA,test_horizon))
cf_WC <- c(rep(NA,150),corf_WC)
plot(c_WC,type = "l",main="Correlation between Wheat and Corn Futures")
lines(cf_WC,type = "l", col = "orange")

c_WS <- c(tail(cor1[1,3,],150),rep(NA,test_horizon))
cf_WS <- c(rep(NA,150),corf_WS)
plot(c_WS,type = "l",main="Correlation between Wheat and Soybeans Futures")
lines(cf_WS,type = "l", col = "orange")

c_CS <- c(tail(cor1[2,3,],150),rep(NA,test_horizon))
cf_CS <- c(rep(NA,150),corf_CS)
plot(c_CS,type = "l",main="Correlation between Corn and Soybeans Futures")
lines(cf_CS,type = "l", col = "orange")

VIII. ARIMAX Models

Based on the correlation analysis above, wheat is most strongly correlated with corn (and vice versa), while soybeans are also most strongly correlated with corn. Each commodity’s ARIMA model is therefore extended with its most-correlated peer as an exogenous regressor (xreg). Note this uses TSA::arima() explicitly rather than forecast::Arima(), since TSA::arima() is the version that supports xreg the way used here.

wheat.corn <-TSA::arima(wheat.ts, , order=c(1,1,3),xreg = corn.ts)
summary(wheat.corn)

Call:
TSA::arima(x = wheat.ts, order = c(1, 1, 3), xreg = corn.ts)

Coefficients:
         ar1      ma1      ma2      ma3    xreg
      0.4740  -0.4003  -0.0352  -0.0768  0.8569
s.e.  0.1071   0.1070   0.0152   0.0129  0.0170

sigma^2 estimated as 0.01275:  log likelihood = 4775.43,  aic = -9540.87

Training set error measures:

              ME RMSE MAE MPE MAPE
Training set NaN  NaN NaN NaN  NaN

checkresiduals(wheat.corn)

    Ljung-Box test

data:  Residuals from ARIMA(1,1,3)
Q* = 11.011, df = 6, p-value = 0.08805

Model df: 4.   Total lags used: 10

acf2(wheat.corn$residuals)

     [,1] [,2] [,3] [,4]  [,5] [,6] [,7] [,8]  [,9] [,10] [,11] [,12] [,13]
ACF     0    0    0    0 -0.02 0.03    0 0.01 -0.02     0     0     0 -0.01
PACF    0    0    0    0 -0.02 0.03    0 0.01 -0.02     0     0     0 -0.01
     [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
ACF   0.02 -0.01 -0.02 -0.02 -0.02 -0.02 -0.02  0.03  0.01 -0.03 -0.02 -0.02
PACF  0.02 -0.01 -0.02 -0.02 -0.02 -0.02 -0.02  0.03  0.01 -0.02 -0.02 -0.02
     [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37]
ACF   0.04 -0.02  0.01 -0.03  0.01  0.00 -0.02 -0.04 -0.01 -0.02  0.02 -0.02
PACF  0.05 -0.02  0.00 -0.04  0.01  0.01 -0.03 -0.04 -0.01 -0.02  0.02 -0.02
     [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48] [,49]
ACF   0.01 -0.03 -0.04 -0.01  0.02  0.02  0.03  0.01 -0.01     0  0.00     0
PACF  0.01 -0.03 -0.05 -0.01  0.01  0.02  0.02  0.01 -0.01     0 -0.01     0
     [,50] [,51] [,52] [,53] [,54] [,55] [,56] [,57] [,58] [,59] [,60] [,61]
ACF   0.02  0.02  0.02 -0.02 -0.02     0 -0.03  0.02  0.00  0.00  0.02 -0.01
PACF  0.01  0.02  0.01 -0.02 -0.01     0 -0.03  0.01 -0.01  0.01  0.02 -0.01
     [,62] [,63] [,64] [,65] [,66] [,67] [,68] [,69] [,70] [,71] [,72] [,73]
ACF  -0.04  0.04 -0.02     0 -0.03  0.01 -0.02 -0.01  0.03 -0.02  0.02 -0.01
PACF -0.04  0.04 -0.03     0 -0.03  0.02 -0.02 -0.01  0.03 -0.03  0.02 -0.02
     [,74] [,75] [,76] [,77] [,78] [,79] [,80] [,81] [,82] [,83] [,84] [,85]
ACF   0.01  0.02 -0.03     0 -0.02 -0.01     0  0.00     0 -0.04 -0.05  0.01
PACF  0.01  0.02 -0.03     0 -0.02 -0.01     0 -0.01     0 -0.03 -0.05  0.00
     [,86] [,87] [,88] [,89] [,90]
ACF  -0.01  0.03 -0.02 -0.03 -0.03
PACF -0.02  0.03 -0.02 -0.04 -0.01

check_normality_residuals(wheat.corn)

Normality Test Summary
Test	Test.Statistic	P.Value	Interpretation
Anderson-Darling Test	181.273	0	Non-Normal
Jarque-Bera Test	481711.226	0	Non-Normal
Lilliefors Test	0.106	0	Non-Normal
Skewness	1.020	N/A	Right-Skewed
Kurtosis	42.890	N/A	Leptokurtic (Heavy Tails)

corn.wheat <-TSA::arima(corn.ts, , order=c(0,1,2),xreg = wheat.ts)
summary(corn.wheat)

Call:
TSA::arima(x = corn.ts, order = c(0, 1, 2), xreg = wheat.ts)

Coefficients:
         ma1      ma2    xreg
      0.0547  -0.0293  0.3329
s.e.  0.0127   0.0129  0.0067

sigma^2 estimated as 0.004994:  log likelihood = 7712.35,  aic = -15418.69

Training set error measures:

              ME RMSE MAE MPE MAPE
Training set NaN  NaN NaN NaN  NaN

checkresiduals(corn.wheat)

    Ljung-Box test

data:  Residuals from ARIMA(0,1,2)
Q* = 16.497, df = 8, p-value = 0.0358

Model df: 2.   Total lags used: 10

acf2(corn.wheat$residuals)

     [,1] [,2] [,3]  [,4]  [,5] [,6] [,7] [,8]  [,9] [,10] [,11] [,12] [,13]
ACF     0    0    0 -0.02 -0.02 0.02 0.03 0.01 -0.02     0  0.01 -0.02  0.01
PACF    0    0    0 -0.02 -0.02 0.02 0.03 0.00 -0.02     0  0.01 -0.02  0.01
     [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
ACF   0.01  0.00 -0.01 -0.01     0  0.01 -0.01 -0.01  0.02 -0.03 -0.04 -0.01
PACF  0.01  0.01 -0.01 -0.01     0  0.01 -0.01 -0.01  0.02 -0.03 -0.04 -0.01
     [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37]
ACF   0.00 -0.02 -0.02 -0.01 -0.01  0.01  0.03     0 -0.01 -0.02  0.06 -0.04
PACF  0.01 -0.02 -0.02 -0.01 -0.01  0.02  0.03     0 -0.01 -0.02  0.06 -0.04
     [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48] [,49]
ACF  -0.01  0.01 -0.01 -0.01  0.04  0.04  0.03  0.01 -0.01  0.04  0.01  0.01
PACF -0.01  0.01 -0.01 -0.01  0.03  0.04  0.03  0.01 -0.01  0.04  0.01  0.00
     [,50] [,51] [,52] [,53] [,54] [,55] [,56] [,57] [,58] [,59] [,60] [,61]
ACF   0.01  0.03 -0.03 -0.02  0.02  0.00 -0.02  0.01  0.02 -0.03 -0.01 -0.01
PACF  0.00  0.03 -0.03 -0.02  0.02  0.01 -0.02  0.01  0.01 -0.02  0.00 -0.01
     [,62] [,63] [,64] [,65] [,66] [,67] [,68] [,69] [,70] [,71] [,72] [,73]
ACF  -0.04  0.02  0.01  0.01  0.00 -0.01  0.01 -0.03  0.01     0  0.00 -0.01
PACF -0.04  0.02  0.00  0.00  0.01 -0.01  0.02 -0.02  0.01     0  0.01 -0.01
     [,74] [,75] [,76] [,77] [,78] [,79] [,80] [,81] [,82] [,83] [,84] [,85]
ACF  -0.02 -0.01 -0.01     0 -0.01  0.01  0.02  0.01     0 -0.01 -0.02  0.01
PACF -0.02 -0.01 -0.02     0 -0.01  0.01  0.02  0.02     0 -0.01 -0.01  0.01
     [,86] [,87] [,88] [,89] [,90]
ACF  -0.01  0.03 -0.03 -0.04 -0.01
PACF -0.02  0.02 -0.03 -0.05 -0.02

soy.corn <-TSA::arima(soy.ts, , order=c(0,1,0),xreg = corn.ts)
summary(soy.corn)

Call:
TSA::arima(x = soy.ts, order = c(0, 1, 0), xreg = corn.ts)

Coefficients:
        xreg
      1.0248
s.e.  0.0229

sigma^2 estimated as 0.02292:  log likelihood = 2938.41,  aic = -5874.81

Training set error measures:

              ME RMSE MAE MPE MAPE
Training set NaN  NaN NaN NaN  NaN

checkresiduals(soy.corn)

    Ljung-Box test

data:  Residuals from ARIMA(0,1,0)
Q* = 21.294, df = 10, p-value = 0.01913

Model df: 0.   Total lags used: 10

acf2(soy.corn$residuals)

      [,1]  [,2]  [,3]  [,4]  [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
ACF  -0.03 -0.01 -0.01 -0.02 -0.04    0 0.01 0.01 0.01  0.01 -0.01 -0.02 -0.02
PACF -0.03 -0.01 -0.01 -0.02 -0.04    0 0.01 0.00 0.01  0.01 -0.01 -0.02 -0.02
     [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
ACF  -0.02  0.01 -0.03 -0.02     0 -0.01 -0.03  0.02  0.01 -0.02 -0.03 -0.01
PACF -0.02  0.01 -0.03 -0.02     0 -0.01 -0.03  0.02  0.01 -0.02 -0.03 -0.01
     [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37]
ACF  -0.03     0  0.01  0.05  0.01     0 -0.03  0.01  0.02     0  0.02  0.01
PACF -0.03     0  0.00  0.04  0.01     0 -0.03  0.01  0.03     0  0.01  0.01
     [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48] [,49]
ACF  -0.01  0.02  0.02  0.01  0.03  0.02  0.03 -0.02 -0.02 -0.04  0.02     0
PACF -0.01  0.01  0.02  0.01  0.03  0.02  0.03 -0.01 -0.02 -0.04  0.02     0
     [,50] [,51] [,52] [,53] [,54] [,55] [,56] [,57] [,58] [,59] [,60] [,61]
ACF  -0.01  0.02 -0.01 -0.01 -0.01  0.03 -0.02 -0.02  0.01 -0.03 -0.01 -0.03
PACF -0.01  0.01 -0.01  0.00  0.00  0.03 -0.01 -0.02  0.00 -0.03 -0.01 -0.03
     [,62] [,63] [,64] [,65] [,66] [,67] [,68] [,69] [,70] [,71] [,72] [,73]
ACF  -0.01  0.02  0.03 -0.04 -0.01  0.01  0.00  0.02  0.00 -0.01 -0.01  0.01
PACF -0.02  0.01  0.03 -0.04 -0.01  0.01  0.01  0.01 -0.01 -0.02 -0.01 -0.01
     [,74] [,75] [,76] [,77] [,78] [,79] [,80] [,81] [,82] [,83] [,84] [,85]
ACF      0 -0.01  0.01  0.00  0.01  0.01 -0.02 -0.02  0.03 -0.02  0.01 -0.01
PACF     0  0.00  0.01 -0.01  0.01  0.00 -0.02 -0.02  0.03 -0.02  0.00 -0.02
     [,86] [,87] [,88] [,89] [,90]
ACF  -0.01  0.00 -0.02     0  0.02
PACF -0.01  0.01 -0.02     0  0.02

IX. Conclusion

From the analysis above, a few key findings emerge:

• No seasonality: neither decompose()/stl(), spectral analysis, nor its bootstrapped confidence intervals turn up significant periodicity in any of the three series — consistent with futures prices tracking macroeconomic conditions rather than a fixed production calendar.
• ARIMA: the best-fitting models (ARIMA(1,1,3) for wheat, ARIMA(0,1,3) for corn, ARIMA(0,1,0) — a random walk — for soybeans) all have residuals that are stationary but not independent or normal, so their forecast intervals should not be trusted.
• GARCH: ARMA(0,0)+GARCH(1,1) fits all three return series well in-sample, with residuals close to independent and normal. Out-of-sample, however, all three forecasts produced negative R² against the actual Oct–Dec 2024 returns — worse than predicting the mean — so these models were not used for further forecasting.
• Multivariate GARCH: the DCC model shows persistent, time-varying correlation across all three commodities, strongest between wheat and corn.
• ARIMAX: adding the most-correlated peer commodity as an exogenous regressor produces statistically significant coefficients, but residuals still fail independence/normality checks, so the same forecasting caveats apply.

Taken together, the value of this analysis is less about producing a deployable price forecast and more about demonstrating the discipline of the model-selection process itself — including the willingness to conclude that a model class (GARCH, here) isn’t reliable enough to forecast with, rather than reporting an unvalidated result.