Wheat price forecasting in R: compare ARIMA with naïve forecasts using rolling-origin backtesting, MAE, RMSE and prediction intervals for monthly World Bank HRW wheat prices.
Can an autoregressive forecasting model really improve on a surprisingly strong benchmark: the most recently observed price? This worked lesson answers the question with out-of-sample predictions, rather than fit statistics or an attractive trend line.
Our CSV contains a World Bank-derived U.S. hard red winter (HRW) Gulf export wheat benchmark, not Pakistan farm-gate prices. The supplied teaching extract was assembled from two public mirrors plus an independently cross-checked October 2026 observation. It has not been compared row-by-row with the current official Excel workbook. The series is in nominal U.S. dollars per metric tonne, and the benchmark specification changed around July 2020. These limitations matter when interpreting structural shifts and historical backtesting. See 01_data/DATA_SOURCE_AND_DICTIONARY.md.
1. Reproduce exactly what you see
Download the full project, unzip it, and launch RStudio with the project root (the folder containing 01_data, 02_R_script, and 03_Quarto_blog) as the working directory. The demonstration uses base R, with no separate modeling package. The Quarto document additionally needs knitr.
The following hidden setup loads each runnable code section directly from the same .R script supplied with the lesson. The article is therefore not a second, independently modified copy of the analysis.
input_file <- "01_data/world_bank_wheat_hrw_monthly_2000_2026.csv"
output_dir <- "04_outputs"
if (!file.exists(input_file)) {
stop("Input CSV not found. Run the script from the AGRON_2026_003 folder.")
}
dir.create(output_dir, showWarnings = FALSE, recursive = TRUE)2. What does the wheat-price series measure?
Each row represents one month. The target is a benchmark export quotation; it is not a futures settlement, farm revenue, a domestic retail price, or an inflation-adjusted measure. Time order is part of the scientific design: the model cannot use a later month’s outcome when predicting an earlier month.
wheat <- read.csv(input_file, stringsAsFactors = FALSE)
if (!identical(names(wheat), c("date", "wheat_price_usd_per_mt"))) {
stop("Unexpected CSV columns; check the supplied Stage 2 dataset.")
}
wheat$date <- as.Date(wheat$date, format = "%Y-%m-%d")
if (anyNA(wheat$date) || anyNA(wheat$wheat_price_usd_per_mt) ||
!is.numeric(wheat$wheat_price_usd_per_mt) ||
any(wheat$wheat_price_usd_per_mt <= 0) || anyDuplicated(wheat$date)) {
stop("Dates and prices must be unique, complete, numeric and positive.")
}
months <- seq.Date(min(wheat$date), max(wheat$date), by = "month")
if (!identical(wheat$date, months) || nrow(wheat) < 240L) {
stop("Expected a chronologically sorted monthly series with at least 240 values.")
}
prices <- wheat$wheat_price_usd_per_mt
price_summary <- data.frame(
first_month = min(wheat$date), last_month = max(wheat$date),
observations = length(prices), min_usd_mt = min(prices),
median_usd_mt = median(prices), mean_usd_mt = mean(prices),
max_usd_mt = max(prices), last_usd_mt = tail(prices, 1)
)
print(price_summary, row.names = FALSE) first_month last_month observations min_usd_mt median_usd_mt mean_usd_mt
2000-01-01 2026-10-01 322 105.11 219.7476 235.6983
max_usd_mt last_usd_mt
522.29 317write.csv(price_summary, file.path(output_dir, "01_price_summary.csv"), row.names = FALSE)The imported dataset has 322 complete monthly observations from January 2000 to October 2026. Its mean price is 235.70 USD/tonne, the observed minimum is 105.11 USD/tonne, and the observed maximum is 522.29 USD/tonne. The last recorded benchmark price is 317.00 USD/tonne. These are observations, not forecasted values. The data come from mixed publication vintages, so they do not reconstruct the historical releases actually available to forecasters at each date.
png(file.path(output_dir, "01_price_history.png"), width = 1700, height = 850, res = 170)
par(mar = c(4.5, 5, 3.4, 1))
plot(wheat$date, prices, type = "l", lwd = 2.2, col = "#19547A",
xlab = "Month", ylab = "Nominal USD per metric tonne",
main = "Monthly U.S. HRW wheat export benchmark, 2000–2026")
abline(v = as.Date("2020-07-01"), col = "#B56B20", lty = 2, lwd = 2)
legend("topleft", "Benchmark specification changes around July 2020",
lty = 2, lwd = 2, col = "#B56B20", bty = "n", cex = 0.83)
box(bty = "l")
dev.off()png
2 
Notice that the range of prices changes substantially through time. A model estimated on calm periods may have difficulty during large market movements. Nothing in this line graph by itself establishes predictability.
3. The experimental question: forecast before seeing the answer
A naïve forecast predicts that the next price will equal the latest available price. It is genuinely useful in financial and commodity series, where the level can move unexpectedly. ARIMA(1,1,0) instead examines the changes in monthly price and allows the most recent change to predict the next one. Here the ARIMA order is fixed before evaluation rather than selected by peeking at outcomes.
The training window grows from December 2016 onward. At every forecast origin, we fit to all observations available up to that month only and predict one and three months ahead. Evaluation uses the same target months (March 2017–October 2026) for both horizons, yielding paired comparisons. This is expanding-window rolling-origin backtesting, not random splitting.
# Each origin uses only prices already observed at that month.
# Fixed ARIMA(1,1,0) is chosen before evaluation, avoiding test-set tuning.
first_origin <- match(as.Date("2016-12-01"), wheat$date)
if (is.na(first_origin)) stop("The required first forecast origin is missing.")
backtest_parts <- lapply(seq.int(first_origin, nrow(wheat) - 1L), function(i) {
training <- prices[seq_len(i)]
fit <- stats::arima(training, order = c(1, 1, 0), method = "ML")
prediction <- stats::predict(fit, n.ahead = 3L)
horizons <- c(1L, 3L)
horizons <- horizons[i + horizons <= nrow(wheat)]
data.frame(
origin_date = wheat$date[i], target_date = wheat$date[i + horizons],
horizon_months = horizons, actual = prices[i + horizons],
naive = rep(tail(training, 1), length(horizons)),
arima = as.numeric(prediction$pred[horizons]),
arima_se = as.numeric(prediction$se[horizons])
)
})
backtest <- do.call(rbind, backtest_parts)
# Identical target months at both horizons make comparisons easier to interpret.
backtest <- subset(backtest, target_date >= as.Date("2017-03-01"))
backtest <- backtest[order(backtest$target_date, backtest$horizon_months), ]
rownames(backtest) <- NULL
if (anyNA(backtest) || any(!is.finite(backtest$arima)) ||
any(backtest$arima_se <= 0) ||
!identical(as.integer(table(backtest$horizon_months)), c(116L, 116L))) {
stop("Unexpected forecast grid or invalid ARIMA predictions.")
}
write.csv(backtest, file.path(output_dir, "02_backtest_forecasts.csv"), row.names = FALSE)
cat("Forecast origins:", length(backtest_parts), " | Paired target months:",
sum(backtest$horizon_months == 1L), "\n")Forecast origins: 118 | Paired target months: 116 The backtest yields 232 predictions: one at each horizon for 116 shared target months. Every scored forecast uses an earlier origin: one month before the target at the short horizon, three months before the target at the longer horizon. We retain the model’s conditional standard error at that origin. This prevents target leakage, though it does not solve the separate issue of later-revised source data.
4. Compare real forecasting errors
For each prediction, the error is the observed price minus the predicted price. Mean absolute error (MAE) averages absolute mistakes in dollars per tonne:
\[ \mathrm{MAE} = \frac{1}{n}\sum_{t=1}^{n}|y_t-\widehat{y}_t|. \]
Root mean squared error (RMSE) gives especially large mistakes more weight:
\[ \mathrm{RMSE} = \sqrt{\frac{1}{n}\sum_{t=1}^{n}(y_t-\widehat{y}_t)^2}. \]
Both should be evaluated on future observations, not the original training sample. Smaller error is better for these scores; neither score is a causal explanation of wheat prices.
accuracy_parts <- lapply(c(1L, 3L), function(h) {
observed <- subset(backtest, horizon_months == h)
rbind(
data.frame(horizon_months = h, model = "Naive", n = nrow(observed),
mae = mean(abs(observed$actual - observed$naive)),
rmse = sqrt(mean((observed$actual - observed$naive)^2))),
data.frame(horizon_months = h, model = "ARIMA(1,1,0)", n = nrow(observed),
mae = mean(abs(observed$actual - observed$arima)),
rmse = sqrt(mean((observed$actual - observed$arima)^2)))
)
})
accuracy <- do.call(rbind, accuracy_parts)
rownames(accuracy) <- NULL
print(accuracy, row.names = FALSE, digits = 4) horizon_months model n mae rmse
1 Naive 116 13.02 19.05
1 ARIMA(1,1,0) 116 13.45 18.94
3 Naive 116 24.59 35.58
3 ARIMA(1,1,0) 116 24.74 36.34write.csv(accuracy, file.path(output_dir, "03_accuracy_by_horizon.csv"), row.names = FALSE)
# Gaussian intervals are conditional approximations, not guaranteed coverage.
z80 <- qnorm(0.90)
z95 <- qnorm(0.975)
backtest$within_80 <- abs(backtest$actual - backtest$arima) <= z80 * backtest$arima_se
backtest$within_95 <- abs(backtest$actual - backtest$arima) <= z95 * backtest$arima_se
coverage <- aggregate(cbind(within_80, within_95) ~ horizon_months,
data = backtest, FUN = mean)
names(coverage)[2:3] <- c("empirical_80", "empirical_95")
coverage$n <- as.integer(table(backtest$horizon_months))
print(coverage, row.names = FALSE, digits = 3) horizon_months empirical_80 empirical_95 n
1 0.836 0.957 116
3 0.853 0.948 116write.csv(coverage, file.path(output_dir, "04_interval_coverage.csv"), row.names = FALSE)
write.csv(backtest, file.path(output_dir, "02_backtest_forecasts.csv"), row.names = FALSE)What did the held-out months actually show?
| Horizon (months) | Model | Test months | MAE | RMSE |
|---|---|---|---|---|
| 1 | Naive | 116 | 13.02 | 19.05 |
| 1 | ARIMA(1,1,0) | 116 | 13.45 | 18.94 |
| 3 | Naive | 116 | 24.59 | 35.58 |
| 3 | ARIMA(1,1,0) | 116 | 24.74 | 36.34 |
Across 116 matched target months, the naïve model has the smaller MAE at both horizons. At one month, its MAE is 13.02, versus 13.45 for ARIMA. ARIMA’s one-month RMSE is marginally lower: 18.94 versus 19.05.
At three months, the naïve method is marginally better on both metrics: MAE 24.59 versus 24.74, and RMSE 35.58 versus 36.34. These differences are small; without a paired uncertainty analysis and an operationally meaningful threshold, they should not be described as a decisive victory. Forecast errors at overlapping origins are dependent, so an ordinary independent-sample significance test would be inappropriate.
one_month <- subset(backtest, horizon_months == 1L &
target_date >= as.Date("2023-01-01"))
png(file.path(output_dir, "02_backtest_comparison.png"), width = 1700, height = 850, res = 170)
par(mar = c(4.5, 5, 3.4, 1))
plot(one_month$target_date, one_month$actual, type = "l", lwd = 2.5,
col = "#202F3A", xlab = "Target month", ylab = "USD per metric tonne",
main = "Real out-of-sample one-month forecasts (recent years)")
lines(one_month$target_date, one_month$naive, col = "#8B9BA4", lwd = 1.6)
lines(one_month$target_date, one_month$arima, col = "#CD7436", lwd = 1.9)
legend("topright", c("Actual", "Naive", "ARIMA(1,1,0)"),
col = c("#202F3A", "#8B9BA4", "#CD7436"),
lwd = c(2.5, 1.6, 1.9), bty = "n")
box(bty = "l")
dev.off()png
2 
mae_table <- xtabs(mae ~ model + horizon_months, data = accuracy)
png(file.path(output_dir, "03_error_metrics.png"), width = 1500, height = 850, res = 170)
par(mar = c(5, 5, 3.5, 1))
barplot(mae_table, beside = TRUE, col = c("#CD7436", "#8B9BA4"),
names.arg = c("1 month", "3 months"),
ylab = "Out-of-sample MAE (USD per metric tonne)",
main = "Prediction errors: smaller is better",
ylim = c(0, max(mae_table) * 1.25), legend.text = TRUE,
args.legend = list(x = "topleft", bty = "n"))
box(bty = "l")
dev.off()png
2 
Should we trust the prediction intervals?
The ARIMA forecast comes with a conditional standard error. We form approximate Gaussian 80% and 95% intervals and ask how often the actual held-out price fell inside them. Their nominal percentages are targets under assumptions, not promises about any single month. Parameter-estimation uncertainty is not fully included in stats::predict.Arima() standard errors, and large shocks can reduce calibration.
Across the same 116 target months per horizon, observed 80% interval coverage is 83.6% at one month (97 / 116) and 85.3% at three months (99 / 116). For the nominal 95% intervals, the corresponding rates are 95.7% and 94.8%.
Observed coverage is near or above the nominal levels in this sample, but this is not proof that the model is correctly specified. Wide intervals can cover many outcomes while having limited decision value. Overlapping forecast origins make coverage outcomes dependent, and this exercise lacks archived original data-release vintages.
5. Fit the full series for the next three months
Only after evaluating historical prediction behavior do we refit the same fixed ARIMA specification to the complete observed series. These forecasts are for the three months after October 2026. They are conditional scenarios—not guaranteed prices or production advice.
final_fit <- stats::arima(prices, order = c(1, 1, 0), method = "ML")
final_prediction <- stats::predict(final_fit, n.ahead = 3L)
future_months <- seq.Date(max(wheat$date), by = "month", length.out = 4L)[-1L]
future <- data.frame(
date = future_months, forecast_usd_mt = as.numeric(final_prediction$pred),
lower_80 = as.numeric(final_prediction$pred - z80 * final_prediction$se),
upper_80 = as.numeric(final_prediction$pred + z80 * final_prediction$se),
lower_95 = as.numeric(final_prediction$pred - z95 * final_prediction$se),
upper_95 = as.numeric(final_prediction$pred + z95 * final_prediction$se)
)
print(future, row.names = FALSE, digits = 5) date forecast_usd_mt lower_80 upper_80 lower_95 upper_95
2026-11-01 310.77 288.63 332.92 276.91 344.64
2026-12-01 309.48 274.76 344.20 256.39 362.58
2027-01-01 309.21 264.80 353.62 241.29 377.13write.csv(future, file.path(output_dir, "05_future_arima_forecast.csv"), row.names = FALSE)The prospective forecast is not an observed price
| Month | Point (USD/tonne) | 80% lower | 80% upper | 95% lower | 95% upper |
|---|---|---|---|---|---|
| Nov 2026 | 310.77 | 288.63 | 332.92 | 276.91 | 344.64 |
| Dec 2026 | 309.48 | 274.76 | 344.20 | 256.39 | 362.58 |
| Jan 2027 | 309.21 | 264.80 | 353.62 | 241.29 | 377.13 |
For November 2026, the model estimates 310.77 USD/tonne. It forecasts 309.21 USD/tonne by January 2027. The point forecast flattens because the effect of the final monthly price change attenuates in this specification. These predictions are not observed prices, and the interval bounds are conditional model-based approximations.
6. Check what the fitted model leaves unexplained
Residuals are the differences between fitted one-step predictions and the corresponding observations. Residual autocorrelation may indicate that the model leaves predictable structure behind. The Ljung–Box test collects serial correlation across the first twelve lags, with one fitted AR parameter accounted for in its degrees of freedom.
A small p-value is evidence against the simplified “no residual serial correlation” null under test assumptions. A large p-value is not proof that the model is accurate, that its intervals are calibrated, or that sudden market shocks are predictable.
residual_values <- as.numeric(stats::residuals(final_fit))
residual_test <- stats::Box.test(residual_values, lag = 12L,
type = "Ljung-Box", fitdf = 1L)
cat("Ljung–Box lag-12 residual check: statistic =",
round(unname(residual_test$statistic), 3),
"; p =", signif(residual_test$p.value, 4), "\n")Ljung–Box lag-12 residual check: statistic = 27.144 ; p = 0.004371 writeLines(c(capture.output(print(final_fit)), "",
capture.output(print(residual_test)), "",
capture.output(sessionInfo())),
file.path(output_dir, "06_model_diagnostics.txt"))
png(file.path(output_dir, "05_residual_diagnostics.png"), width = 1600, height = 760, res = 170)
par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
plot(wheat$date, residual_values, type = "l", col = "#19547A",
xlab = "Month", ylab = "ARIMA residual (USD/mt)",
main = "Full-series model residuals")
abline(h = 0, col = "#CD7436", lty = 2)
stats::acf(residual_values, main = "Residual autocorrelation")
dev.off()png
2 The fitted AR coefficient is 0.2076 (with reported standard error 0.0547). This describes dependence between consecutive price changes under the model; it is neither a causal parameter nor proof of predictive superiority. The Ljung–Box test at 12 lags returns a p-value of 0.004371. For the uploaded execution it is below 0.01, indicating remaining residual serial dependence. Consequently, the fitted model should not be called fully adequate merely because its interval coverage is close to nominal. Large residual spikes also suggest periods of changing volatility and potential structural instability.
recent <- subset(wheat, date >= as.Date("2023-01-01"))
png(file.path(output_dir, "04_three_month_forecast.png"), width = 1600, height = 840, res = 170)
par(mar = c(4.5, 5, 3.4, 1))
plot(c(recent$date, future$date),
c(recent$wheat_price_usd_per_mt, future$forecast_usd_mt),
type = "n", xlab = "Month", ylab = "USD per metric tonne",
ylim = range(c(recent$wheat_price_usd_per_mt, future$lower_95,
future$upper_95)),
main = "Next three months: ARIMA forecast and conditional intervals")
polygon(c(tail(recent$date, 1), future$date, rev(future$date)),
c(tail(prices, 1), future$lower_95, rev(future$upper_95)),
col = "#DCE7EF", border = NA)
polygon(c(tail(recent$date, 1), future$date, rev(future$date)),
c(tail(prices, 1), future$lower_80, rev(future$upper_80)),
col = "#A7C5DB", border = NA)
lines(recent$date, recent$wheat_price_usd_per_mt,
col = "#202F3A", lwd = 2.3)
lines(c(tail(recent$date, 1), future$date),
c(tail(prices, 1), future$forecast_usd_mt),
col = "#CD7436", lwd = 2.2, type = "b", pch = 16)
legend("topleft", c("Observed", "ARIMA point forecast", "80% interval", "95% interval"),
col = c("#202F3A", "#CD7436", "#A7C5DB", "#DCE7EF"),
lwd = c(2.3, 2.2, 7, 7), bty = "n", cex = 0.8)
box(bty = "l")
dev.off()png
2 
cat("Wrote", length(list.files(output_dir)), "files to", output_dir, "\n")Wrote 13 files to 04_outputs 7. Interpreting the evidence without overclaiming
Observed conclusion: In the supplied 2017–2026 matched-target backtest, ARIMA(1,1,0) does not establish a dependable improvement over carrying the last observed price forward. The naïve method wins on MAE at both horizons and on RMSE at three months; ARIMA’s one-month RMSE advantage is small. The model-based 80% and 95% interval coverage rates are descriptively near nominal, yet the Ljung–Box check detects leftover serial correlation. These findings show why predictive accuracy, calibration and residual adequacy must be assessed separately.
The narrow comparison is between two prespecified methods on one teaching dataset and a single backtest design. It does not rank every ARIMA variant, forecast origin window or agricultural commodity.
Key limitations: (1) the dataset merges published sources/vintages rather than replaying actual historical data releases; (2) the HRW benchmark definition changes around July 2020; (3) monthly quotes are nominal and do not represent farm prices; (4) no external covariates such as weather, crop supply, policy or exchange rates are included; (5) substantial market shocks may violate model assumptions. These models describe associations in a time series, not causes.
Adapt this to a thesis: replace the CSV with a consistently defined official monthly series; document revisions and units; register the proposed models; reserve genuinely later dates for evaluation; compare multiple horizons; examine residuals and empirical interval coverage; and explain when the model fails. For real operational use, backtest against the data releases that were actually available at each forecast origin, if historical vintages can be obtained.
Data, reproducibility and citations
- Primary source: World Bank Commodity Markets (Pink Sheet), series
WHEAT_US_HRW, monthly nominal USD/metric tonne. - Source catalog and terms: World Bank Commodity Prices catalog.
- Dataset acquisition: mirror A, mirror B, and an independently cross-checked October 2026 quote. Read the bundled documentation for details and precise blob fingerprints.
- R statistical reference: R
stats::predict.Arima()documentation. - Method: Fixed ARIMA(1,1,0) using base R
stats::arima(); rolling origins; paired target months; MAE, RMSE, and out-of-sample interval coverage.
Stage 4 evidence and execution note: All outcome-specific claims were checked against user-uploaded real R output files dated 7 November 2026; accuracy and empirical coverage were independently recalculated from the 232-row forecast table. The diagnostic text reports R 4.5.1, and the supplied original HTML identifies Quarto 1.8.25. This revised Quarto source itself has not been freshly rendered in the current environment (R and Quarto are unavailable here). Its inline values and tables will be recomputed when the user renders it locally. The already-supplied Stage 3 standalone HTML may not reflect these newer result-specific interpretations until re-rendered.
Comments
Post a Comment