Thailand Agriculture Export Data

Library

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.4.3
library(ggfortify)
## Warning: package 'ggfortify' was built under R version 4.4.3
library(fpp2)
## Warning: package 'fpp2' was built under R version 4.4.3
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## Registered S3 methods overwritten by 'forecast':
##   method                 from     
##   autoplot.Arima         ggfortify
##   autoplot.acf           ggfortify
##   autoplot.ar            ggfortify
##   autoplot.bats          ggfortify
##   autoplot.decomposed.ts ggfortify
##   autoplot.ets           ggfortify
##   autoplot.forecast      ggfortify
##   autoplot.stl           ggfortify
##   autoplot.ts            ggfortify
##   fitted.ar              ggfortify
##   fortify.ts             ggfortify
##   residuals.ar           ggfortify
## ── Attaching packages ────────────────────────────────────────────── fpp2 2.5 ──
## ✔ forecast  8.23.0     ✔ expsmooth 2.3   
## ✔ fma       2.5
## Warning: package 'forecast' was built under R version 4.4.3
## Warning: package 'fma' was built under R version 4.4.3
## Warning: package 'expsmooth' was built under R version 4.4.3
## 
library(forecast)
library(TSstudio)
## Warning: package 'TSstudio' was built under R version 4.4.3
library(urca)
## Warning: package 'urca' was built under R version 4.4.3
library(tseries)
## Warning: package 'tseries' was built under R version 4.4.3

Data

data <- read.csv('AgriData.csv')

export_data <- ts(data$Export,frequency = 12,start = c(2013,1))

export_data
##         Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep    Oct
## 2013 109709 102289 113264  98807 103899  96357 104308 102112 102517 113344
## 2014 108612 107994 111678  99035 110895 108042 111282 106996 107398 118103
## 2015  91847  89265 105526  94948 106776 100471 103255 101373 102246 110696
## 2016  94939  98230 113115  95754 104244 100579  89340  99964 100228  99474
## 2017 101704 110668 127149 107441 122187 113967 111223 119655 108031 105331
## 2018 108601 102450 108701 101148 109728 109028 105886 117532 105272 115007
## 2019 109730 106398 116954 117303 118257 104236 109850 115373 101271 110027
## 2020  96777 100944 117541 127799 125817 100201 101117 104139 106319 102294
## 2021  98371 105125 119388 128650 136497 130283 126405 132736 126691 128570
## 2022 123726 131230 146181 155616 187118 178900 163144 157075 141754 141481
## 2023 126686 136754 164767 172943 160342 162490 142872 141951 150505 145313
##         Nov    Dec
## 2013 107799 113934
## 2014 103104 114799
## 2015  98434 108092
## 2016 109125 113014
## 2017 124026 112869
## 2018 110090 103718
## 2019 106619 101914
## 2020 104857 102237
## 2021 131393 138903
## 2022 152378 137815
## 2023 147212 137815

Time Series Decomposition

autoplot(export_data) + labs(title = "Thai agricultural export value", y = "Million baht", x = NULL)

ggseasonplot(export_data,continuous = TRUE,polar = TRUE)
## Warning in fortify(data, ...): Arguments in `...` must be used.
## ✖ Problematic argument:
## • na.rm = TRUE
## ℹ Did you misspell an argument name?

ggsubseriesplot(export_data)
## Warning in fortify(data, ...): Arguments in `...` must be used.
## ✖ Problematic argument:
## • na.rm = TRUE
## ℹ Did you misspell an argument name?

ts_seasonal(export_data,type = "all")

White noise?

From Ljung-Box we can see that the data isn’t white-noise so it has something to forecast.

Box.test(export_data, lag = 24, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  export_data
## X-squared = 843.51, df = 24, p-value < 2.2e-16

Stationarity (For ARIMA model)

After one diff you can see that the time series is stationary so we can show that optimal params for ARIMA is d=1 and D=0

#Unit Root Test

export_data %>% diff() %>% ur.kpss() %>% summary() #KPSS
## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 4 lags. 
## 
## Value of test-statistic is: 0.0469 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739
export_data %>% diff() %>% adf.test() #ADF
## Warning in adf.test(.): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  .
## Dickey-Fuller = -6.1381, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary
export_data %>% diff() %>% ggtsdisplay()

nsdiffs(export_data)
## [1] 0
BoxCox.lambda(export_data)
## [1] -0.9999242

Train Test Split

We use data from 2013- 9/2022 to forecast to avoid unexpected growth in trend.

n <- length(export_data)              
n_train <- floor(0.8 * n)             

train <- window(export_data, end = c(2022, 9))
test  <- window(export_data, start = c(2022, 10))
h <- length(test)                     

length(train); length(test)                           
## [1] 117
## [1] 15

Model Selection

In this step we set Naive Forecasting Method to the baseline.

fc_naive  <- naive(train, h = h)
fc_snaive <- snaive(train, h = h)

# ETS 
fit_ets <- ets(train)
fc_ets  <- forecast(fit_ets, h = h)

# Holt-Winters, additive and multiplicative seasonality
fc_hw_add  <- hw(train, h = h, seasonal = "additive")
fc_hw_mult <- hw(train, h = h, seasonal = "multiplicative")

# ARIMA and SARIMA,
fit_arima  <- auto.arima(train, seasonal = FALSE, stepwise = FALSE, approximation = FALSE)
fit_sarima <- auto.arima(train, seasonal = TRUE,  stepwise = FALSE, approximation = FALSE)
fc_arima   <- forecast(fit_arima,  h = h)
fc_sarima  <- forecast(fit_sarima, h = h)

Residual Checks

Since SARIMA is the only model that error is white-noise so we will use it for point and interval forecast.

lb <- function(res, fitdf = 0) {
  r <- Box.test(na.omit(res), lag = 24, type = "Ljung-Box", fitdf = fitdf)
  c(Q = unname(r$statistic), p = r$p.value)
}

resids <- list(
  Naive   = residuals(fc_naive),
  SNaive  = residuals(fc_snaive),
  ETS     = residuals(fit_ets),
  HW_add  = residuals(fc_hw_add),
  HW_mult = residuals(fc_hw_mult),
  ARIMA   = residuals(fit_arima),
  SARIMA  = residuals(fit_sarima)
)
fitdf <- c(Naive = 0, SNaive = 0, ETS = 0, HW_add = 0, HW_mult = 0,
           ARIMA = length(coef(fit_arima)), SARIMA = length(coef(fit_sarima)))

round(t(mapply(lb, resids, fitdf)), 4)
##                Q      p
## Naive   139.8440 0.0000
## SNaive  280.8548 0.0000
## ETS      61.1297 0.0000
## HW_add   58.2162 0.0001
## HW_mult  55.1995 0.0003
## ARIMA    49.0298 0.0002
## SARIMA   18.5411 0.4866
checkresiduals(fit_ets, lag = 24)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(M,N,A)
## Q* = 61.13, df = 24, p-value = 4.417e-05
## 
## Model df: 0.   Total lags used: 24
checkresiduals(fit_arima,lag = 24)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,5)
## Q* = 49.03, df = 19, p-value = 0.0001818
## 
## Model df: 5.   Total lags used: 24
checkresiduals(fit_sarima, lag = 24)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,4)(1,0,0)[12]
## Q* = 18.541, df = 19, p-value = 0.4866
## 
## Model df: 5.   Total lags used: 24

Model Comparison

metrics <- function(pred, actual, train) {
  e <- as.numeric(actual) - as.numeric(pred)
  c(RMSE = sqrt(mean(e^2)),
    MAE  = mean(abs(e)),
    MAPE = mean(abs(e / as.numeric(actual))) * 100)
}

preds <- list(Naive = fc_naive$mean, SNaive = fc_snaive$mean,
              ETS = fc_ets$mean, HW_add = fc_hw_add$mean, HW_mult = fc_hw_mult$mean,
              ARIMA = fc_arima$mean, SARIMA = fc_sarima$mean)

res <- t(sapply(preds, metrics, actual = test, train = train))
round(res[order(res[, "RMSE"]), ], 2)
##             RMSE      MAE  MAPE
## ETS     10277.31  8716.16  5.88
## HW_add  10499.60  8981.74  6.05
## Naive   13648.42 10096.80  6.51
## SNaive  15335.04 13357.60  8.79
## ARIMA   15852.30 13566.75  9.57
## SARIMA  19050.05 16412.12 11.38
## HW_mult 24985.15 20880.82 14.16
autoplot(window(export_data, start = 2020))+
  autolayer(fc_snaive$mean,series = "SNaive")+
  autolayer(fc_ets$mean,series = "ETS")+
  autolayer(fc_hw_add$mean,series = "HW additive")+
  autolayer(fc_sarima$mean,series = "SARIMA")+
  autolayer(test,series = "Actual")

## Point and Interval Forecast

# Check that the intervals behave on the holdout
fc <- forecast(fit_sarima, h = h, level = c(80, 95))
mean(test >= fc$lower[, 2] & test <= fc$upper[, 2])   
## [1] 0.8666667
mean(test >= fc$lower[, 1] & test <= fc$upper[, 1])   
## [1] 0.4666667
# Final model on all data, and the 12-month forecast
fit_final <- auto.arima(export_data, stepwise = FALSE, approximation = FALSE)
checkresiduals(fit_final, lag = 24)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(4,1,0)(1,0,0)[12]
## Q* = 13.221, df = 19, p-value = 0.8271
## 
## Model df: 5.   Total lags used: 24
fc_final <- forecast(fit_final, h = 12, level = c(80, 95))
fc_final
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jan 2024       129003.6 118799.1 139208.1 113397.2 144610.0
## Feb 2024       137444.6 124855.2 150034.0 118190.7 156698.5
## Mar 2024       153462.7 139427.0 167498.4 131996.9 174928.5
## Apr 2024       157806.1 142657.2 172955.0 134637.9 180974.3
## May 2024       151270.6 135708.5 166832.8 127470.4 175070.9
## Jun 2024       151368.7 134958.0 167779.4 126270.7 176466.7
## Jul 2024       140835.8 123509.1 158162.4 114337.0 167334.6
## Aug 2024       140534.5 122332.5 158736.5 112697.0 168372.0
## Sep 2024       145148.3 126018.0 164278.5 115891.0 174405.5
## Oct 2024       142624.4 122777.3 162471.4 112271.0 172977.8
## Nov 2024       143546.8 123006.5 164087.2 112133.0 174960.6
## Dec 2024       138435.9 117208.1 159663.7 105970.8 170901.0
autoplot(fc_final)

## Conclusion

Thai agricultural export value is a non-stationary, seasonal series that needs one difference (d = 1) and no variance transformation. Of the seven models compared, SARIMA(0,1,4)(1,0,0)[12] is the final choice because it’s the only one with valid residuals, so it’s the only one whose prediction intervals you can defend.

what we find in this project:

1.ETS was best (RMSE 10,277, MAPE 5.88%), followed closely by Holt-Winters additive (10,500). Naive scored 13,648, and SARIMA scored 19,050, worse than both Naive and ETS.

2.Residual diagnostics: only SARIMA passed the Ljung-Box test (p = 0.49). ETS and Holt-Winters failed because they miss the lag-4 and lag-12 dependence.

3.Final model: SARIMA, chosen for statistical validity, with the forecast and its 80% and 95% intervals taken from that one model.