This example fitted the STCCGEV model for regions Dufferin and Wellington with covariates cdd, frost_days, rx1day, tg_mean, and txgt_25 and predicted crop yields and compared them with actual data.
``` r
bsts_Dufferin <- fit_bsts(yy_train[,1], zz_train[,1,], lags = 2, MCMC.iter = 10)
#> =-=-=-=-= Iteration 0 Sun Sep 27 09:56:19 2026 =-=-=-=-=
#> =-=-=-=-= Iteration 1 Sun Sep 27 09:56:19 2026 =-=-=-=-=
#> =-=-=-=-= Iteration 2 Sun Sep 27 09:56:19 2026 =-=-=-=-=
#> =-=-=-=-= Iteration 3 Sun Sep 27 09:56:19 2026 =-=-=-=-=
#> =-=-=-=-= Iteration 4 Sun Sep 27 09:56:19 2026 =-=-=-=-=
#> =-=-=-=-= Iteration 5 Sun Sep 27 09:56:19 2026 =-=-=-=-=
#> =-=-=-=-= Iteration 6 Sun Sep 27 09:56:19 2026 =-=-=-=-=
#> =-=-=-=-= Iteration 7 Sun Sep 27 09:56:19 2026 =-=-=-=-=
#> =-=-=-=-= Iteration 8 Sun Sep 27 09:56:19 2026 =-=-=-=-=
#> =-=-=-=-= Iteration 9 Sun Sep 27 09:56:19 2026 =-=-=-=-=
bsts_Wellington <- fit_bsts(yy_train[,2], zz_train[,2,], lags = 2, MCMC.iter = 10)
#> =-=-=-=-= Iteration 0 Sun Sep 27 09:56:19 2026 =-=-=-=-=
#> =-=-=-=-= Iteration 1 Sun Sep 27 09:56:19 2026 =-=-=-=-=
#> =-=-=-=-= Iteration 2 Sun Sep 27 09:56:19 2026 =-=-=-=-=
#> =-=-=-=-= Iteration 3 Sun Sep 27 09:56:19 2026 =-=-=-=-=
#> =-=-=-=-= Iteration 4 Sun Sep 27 09:56:19 2026 =-=-=-=-=
#> =-=-=-=-= Iteration 5 Sun Sep 27 09:56:19 2026 =-=-=-=-=
#> =-=-=-=-= Iteration 6 Sun Sep 27 09:56:19 2026 =-=-=-=-=
#> =-=-=-=-= Iteration 7 Sun Sep 27 09:56:19 2026 =-=-=-=-=
#> =-=-=-=-= Iteration 8 Sun Sep 27 09:56:19 2026 =-=-=-=-=
#> =-=-=-=-= Iteration 9 Sun Sep 27 09:56:19 2026 =-=-=-=-=
list_bsts_sample <- list(bsts_Dufferin, bsts_Wellington)
Gaussianforecasts_G <- simulation_generalized(nsim = 10,
n_train = n_train,
n_test = n_test,
copula = "Gaussian",
init_params = init_params_full_G,
fn = log_likelihood_Generalized,
U_train = uu,
Z_train = zz_train,
X = xx_train,
Y_test = yy_test,
BSTS_list = list_bsts_sample)
Dufferin_Gaussian_plot<- plot_forecast(forecast = Gaussianforecasts_G[[3]][,,1],
data_train = yy_train[,1],
data_test = yy_test[,1],
time = time_all,
quant_high = 0.95,
quant_low = 0.05,
observed_col = "#e23345",
forecast_col = "#CF9FFF",
title = "Dufferin - Gaussian copula forecast")
Wellington_Gaussian_plot<- plot_forecast(forecast = Gaussianforecasts_G[[3]][,,2],
data_train = yy_train[,2],
data_test = yy_test[,2],
time = time_all,
quant_high = 0.95,
quant_low = 0.05,
observed_col = "#6195c4",
forecast_col = "#CF9FFF",
title = "Wellington - Gaussian copula forecast")
print(Dufferin_Gaussian_plot)