This example demonstrates the fitting of a spatio-temporal model to predict crop yields in Chatham-Kent, Ontario, and compares the predictions to real data.
``` r
bsts_1 <- fit_bsts(y_train[[medoid_names[1]]], z_train[[medoid_names[1]]], lags = 2, MCMC.iter = 10) # Chatham-Kent
#> =-=-=-=-= Iteration 0 Thu Sep 4 16:18:40 2025 =-=-=-=-=
#> =-=-=-=-= Iteration 1 Thu Sep 4 16:18:40 2025 =-=-=-=-=
#> =-=-=-=-= Iteration 2 Thu Sep 4 16:18:40 2025 =-=-=-=-=
#> =-=-=-=-= Iteration 3 Thu Sep 4 16:18:40 2025 =-=-=-=-=
#> =-=-=-=-= Iteration 4 Thu Sep 4 16:18:40 2025 =-=-=-=-=
#> =-=-=-=-= Iteration 5 Thu Sep 4 16:18:40 2025 =-=-=-=-=
#> =-=-=-=-= Iteration 6 Thu Sep 4 16:18:40 2025 =-=-=-=-=
#> =-=-=-=-= Iteration 7 Thu Sep 4 16:18:40 2025 =-=-=-=-=
#> =-=-=-=-= Iteration 8 Thu Sep 4 16:18:40 2025 =-=-=-=-=
#> =-=-=-=-= Iteration 9 Thu Sep 4 16:18:40 2025 =-=-=-=-=
bsts_2 <- fit_bsts(y_train[[medoid_names[2]]], z_train[[medoid_names[2]]], lags = 2, MCMC.iter = 10) # Wellington
#> =-=-=-=-= Iteration 0 Thu Sep 4 16:18:40 2025 =-=-=-=-=
#> =-=-=-=-= Iteration 1 Thu Sep 4 16:18:40 2025 =-=-=-=-=
#> =-=-=-=-= Iteration 2 Thu Sep 4 16:18:40 2025 =-=-=-=-=
#> =-=-=-=-= Iteration 3 Thu Sep 4 16:18:40 2025 =-=-=-=-=
#> =-=-=-=-= Iteration 4 Thu Sep 4 16:18:40 2025 =-=-=-=-=
#> =-=-=-=-= Iteration 5 Thu Sep 4 16:18:40 2025 =-=-=-=-=
#> =-=-=-=-= Iteration 6 Thu Sep 4 16:18:40 2025 =-=-=-=-=
#> =-=-=-=-= Iteration 7 Thu Sep 4 16:18:40 2025 =-=-=-=-=
#> =-=-=-=-= Iteration 8 Thu Sep 4 16:18:40 2025 =-=-=-=-=
#> =-=-=-=-= Iteration 9 Thu Sep 4 16:18:40 2025 =-=-=-=-=
bsts_3 <- fit_bsts(y_train[[medoid_names[3]]], z_train[[medoid_names[3]]], lags = 1, MCMC.iter = 10) # Lambton
#> =-=-=-=-= Iteration 0 Thu Sep 4 16:18:40 2025 =-=-=-=-=
#> =-=-=-=-= Iteration 1 Thu Sep 4 16:18:40 2025 =-=-=-=-=
#> =-=-=-=-= Iteration 2 Thu Sep 4 16:18:40 2025 =-=-=-=-=
#> =-=-=-=-= Iteration 3 Thu Sep 4 16:18:40 2025 =-=-=-=-=
#> =-=-=-=-= Iteration 4 Thu Sep 4 16:18:40 2025 =-=-=-=-=
#> =-=-=-=-= Iteration 5 Thu Sep 4 16:18:40 2025 =-=-=-=-=
#> =-=-=-=-= Iteration 6 Thu Sep 4 16:18:40 2025 =-=-=-=-=
#> =-=-=-=-= Iteration 7 Thu Sep 4 16:18:40 2025 =-=-=-=-=
#> =-=-=-=-= Iteration 8 Thu Sep 4 16:18:40 2025 =-=-=-=-=
#> =-=-=-=-= Iteration 9 Thu Sep 4 16:18:40 2025 =-=-=-=-=
list_bsts <- list(bsts_1, bsts_2, bsts_3)
GaussianForecast3 <- simul_fun_noGEV_3d(nsim = 10,
n_train = n_train,
n_test = n_test,
copula = copula_list[1],
init_params = init_params_full,
fn = log_likelihood_noGEV_3d,
u1 = u[[1]],
u2 = u[[2]],
u3 = u[[3]],
z1_train = z_train[[1]],
z2_train = z_train[[2]],
z3_train = z_train[[3]],
z1_test = z_test[[1]],
z2_test = z_test[[2]],
z3_test = z_test[[3]],
X_t = (z_train[[1]] + z_train[[2]] + z_train[[3]])/3,
y1_test = y_test[[1]],
y2_test = y_test[[2]],
y3_test = y_test[[3]],
BSTS_1 = bsts_1,
BSTS_2 = bsts_2,
BSTS_3 = bsts_3)
ChathamKent_copula_plot <- plot_forecast(forecast = GaussianForecast3[[2]],
data_train = y_train$`Chatham-Kent`,
data_test = y_test$`Chatham-Kent`,
time = time,
quant_high = 0.95,
quant_low = 0.05,
observed_col = "#96BA66",
forecast_col = "#CF9FFF",
title = "Chatham-Kent - Gaussian copula forecast")
print(ChathamKent_copula_plot)