---
title: "Applying Spatio-Temporal Model to Crop Yield Data"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Introduction to STCYP}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r setup}
library(STCYP)
```
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
bsts_2 <- fit_bsts(y_train[[medoid_names[2]]], z_train[[medoid_names[2]]], lags = 2, MCMC.iter = 10) # Wellington
bsts_3 <- fit_bsts(y_train[[medoid_names[3]]], z_train[[medoid_names[3]]], lags = 1, MCMC.iter = 10) # Lambton

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)                              
```