---
title: "GDP Analysis Using Structural Time Series Models: The Multiplicative Case"
author: "Kätlin Kippar"
date: "`r Sys.Date()`"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


***

#### Packages:
```{r message=FALSE, warning=FALSE}
library(ggplot2)
library(autostsm)
library(forecast)
library(dplyr)
library(lubridate)
library(mFilter)
library(lmtest)
library(vars)
```



***

#### Data:
```{r message=FALSE, warning=FALSE}
estonia <- readr::read_delim("CLVMNACSCAB1GQEE.csv")
finland <- readr::read_delim("CLVMNACSCAB1GQFI.csv")
germany <- readr::read_delim("CLVMNACSCAB1GQDE.csv")
latvia <- readr::read_delim("CLVMNACSCAB1GQLV.csv")
lithuania <- readr::read_delim("CLVMNACSCAB1GQLT.csv")
sweden <- readr::read_delim("CLVMNACSCAB1GQSE.csv")
```


```{r}
colnames(estonia) <- c("date", "gdp")

colnames(latvia) <- c("date", "gdp")

colnames(lithuania) <- c("date", "gdp")

colnames(sweden) <- c("date", "gdp")
sweden <- sweden |>  filter(date > as.Date('1994-10-01') )

colnames(finland) <- c("date", "gdp")
finland <- finland |> filter(date > as.Date('1994-10-01') )

colnames(germany) <- c("date", "gdp")
germany <- germany |> filter(date > as.Date('1994-10-01') )
```



***

#### Models:
#### All models:
```{r}
model_fun <- function(df, multp){
  ## Empty data frames to store all results:
  stsm_trig_df <- data.frame()
  stsm_fitted_trig_df <- data.frame()
  stsm_arma_df <- data.frame()
  stsm_fitted_arma_df <- data.frame()
  
  ## Model estimation:
  stsm_trig <- stsm_estimate(df, freq=4, seasons=F,det_drift=F, det_cycle = F, det_trend = F, decomp="trend-cycle", cycle = "trig", multiplicative = multp)
  stsm_arma <- stsm_estimate(df, freq=4, seasons=F,det_drift=F, det_cycle = F, det_trend = F, decomp="trend-cycle", cycle = "arma", multiplicative = multp)
  
  ## Fitting the model:
  stsm_fitted_trig <- stsm_forecast(model=stsm_trig, y=df, n.ahead = 0) 
  stsm_fitted_arma <- stsm_forecast(model=stsm_arma, y=df, n.ahead = 0) 
  
  
  ## Combining results into a single data frame
  stsm_trig_df <- rbind(stsm_trig_df, stsm_trig)
  stsm_fitted_trig_df <- rbind(stsm_fitted_trig_df, stsm_fitted_trig)
  stsm_arma_df <- rbind(stsm_arma_df, stsm_arma)
  stsm_fitted_arma_df <- rbind(stsm_fitted_arma_df, stsm_fitted_arma)
    
  
  return(list(stsm_trig_df=stsm_trig_df, stsm_fitted_trig_df=stsm_fitted_trig_df, stsm_arma_df=stsm_arma_df, stsm_fitted_arma_df=stsm_fitted_arma_df))
}

```



```{r}
models_fun <- function(df_list, multp){
  ## Empty data frames to store all results
  model_form_trig <- data.frame()
  model_trig <- data.frame()
  model_form_arma <- data.frame()
  model_arma <- data.frame()
  
  for (j in 1:length(df_list)){ ## Iterate through the list of GDP data frames for each country
    ## Estimating models:
    m_res <- model_fun(get(df_list[j]), multp)
    
    ## Separating the ARMA and trig models:
    stsm_form_trig <- m_res$stsm_trig_df
    stsm_trig <- m_res$stsm_fitted_trig_df
    stsm_form_arma <- m_res$stsm_arma_df
    stsm_arma <- m_res$stsm_fitted_arma_df
    
    ## Including a variable containing the country name:
    stsm_form_trig <- stsm_form_trig |> mutate(country = as.character(df_list[j]) ) 
    stsm_trig <- stsm_trig |> mutate(country = as.character(df_list[j]) ) 
    stsm_form_arma <- stsm_form_arma |> mutate(country = as.character(df_list[j]) ) 
    stsm_arma <- stsm_arma |> mutate(country = as.character(df_list[j]) ) 
 
    ## Combining results into a single data frame
    model_form_trig <- rbind(model_form_trig, stsm_form_trig)
    model_trig <- rbind(model_trig, stsm_trig)
    model_form_arma <- rbind(model_form_arma, stsm_form_arma)
    model_arma <- rbind(model_arma, stsm_arma)
    
  }
    return(list(model_form_trig=model_form_trig,model_trig=model_trig, model_form_arma=model_form_arma,model_arma=model_arma)) 
}

```


```{r warning=FALSE}
all_models <- models_fun(c("estonia","latvia", "lithuania", "sweden", "finland", "germany"), multp = TRUE) 

trig_models <- all_models$model_trig ## Fitted trigonometric models
trig_model_forms <- all_models$model_form_trig  ## Trigonometric model formulas with estimated coefficients

arma_models <- all_models$model_arma
arma_model_forms <- all_models$model_form_arma
```




***


***



#### Relations Between the Business Cycles 
##### Granger Causality Test

##### Functions:
```{r}
granger_test_fun <- function(CN, type){ ## CN - county name
  model_name <- paste0("model_", type)
  model_df <- all_models[[model_name]]
  model_df_estonia <- model_df |> filter(country=="estonia")
  model_df_CN <-  model_df |> filter(country==CN)
  
  ## Creating the cycle data frame required for variable selection and the Granger causality test:
  cycle_df <- cbind(estonia=(model_df_estonia$cycle),CN=(model_df_CN$cycle) ) |> as.data.frame()
  
  ## Selecting the lag order:
  var_sel <- VARselect(cycle_df, lag.max = 10, type = c("const"))
  
  ## Testing whether the business cycle values of **Country** are causing Estonia’s business cycle values
  gt <- grangertest(estonia~CN, order = var_sel$selection[[1]], data=cycle_df) 
  
  ## Testing whether Estonia’s business cycle values are causing the business cycle values of **Country**
  gt2 <- grangertest(CN~estonia, order = var_sel$selection[[1]], data=cycle_df) 
  
  return(list(order=var_sel$selection[[1]],gt,gt2))
  
}



granger_test_fun_hp <- function(CN){
  hp_estonia <- hpfilter(log(estonia$gdp),freq=1600,type=c("lambda"), drift = T)
  CN_df <- get(CN)
  hp_CN <- hpfilter(log(CN_df$gdp),freq=1600,type=c("lambda"), drift = T)
  
  hp_cycle <-  cbind(estonia=(hp_estonia$cycle),CN=(hp_CN$cycle ) ) |> as.data.frame()
  
  var_sel <- VARselect(hp_cycle, lag.max = 10, type = c("const"))
  gt <- grangertest(estonia~CN, order = var_sel$selection[[1]], data=hp_cycle) 
  gt2 <- grangertest(CN~estonia, order = var_sel$selection[[1]], data=hp_cycle) 
  
  return(list(order=var_sel$selection[[1]],gt,gt2))
  
}

## This function combines the Granger test results for all countries into a single table
fun_all_results_gt <- function(countries, model_type){
    table <- data.frame()
    for (j in 1:length(countries)){
      for (k in 1:length(model_type)){
        res <- granger_test_fun(countries[j], model_type[k])
        order <- res$order
        row <- c(countries[j], model_type[k], order, res[[2]]$`Pr(>F)`[2],  res[[3]]$`Pr(>F)`[2])
        table <- rbind(table, row)
      }
    
      res_hp <-  granger_test_fun_hp(countries[j])
      order_hp <- res_hp$order
      row_hp  <-  c(countries[j], "hp" , order_hp, res_hp[[2]]$`Pr(>F)`[2],  res_hp[[3]]$`Pr(>F)`[2])
      table <- rbind(table, row_hp)
      
      
    }
    colnames(table) <- c("county", "model_type", "order", "p-value_country_causes_estonia", "p-value_estonia_causes_country")

  return(table)
}

```



```{r}
granger_test_res <- fun_all_results_gt(c("latvia", "lithuania", "sweden", "finland", "germany"), c("trig", "arma"))
granger_test_res$`p-value_country_causes_estonia` <- round(as.numeric(granger_test_res$`p-value_country_causes_estonia`),4)
granger_test_res$`p-value_estonia_causes_country` <- round(as.numeric(granger_test_res$`p-value_estonia_causes_country`),4)
granger_test_res 
```


```{r}
# writexl::write_xlsx(granger_test_res, "mltp_granger_test_res.xlsx")
```




***


***




#### Forecasting
##### RMSE and MAPE

##### Function for finding predictions:

The purpose of this function is to calculate forecasts for a given country. It generates predictions for the next six periods, each time expanding the input data by one additional observation. The first forecast is based on data from 1995-01-01 to 2020-10-01.

Note: The first prediction date is computed by subtracting five quarters from the end of the dataset. Since the training set uses < n_date (excluding the observation at n_date), this adjustment ensures that the final 6-step-ahead forecast lands exactly on 2024-10-01, the last available GDP value.
```{r}
prediction_function <- function(country_df, multp){
  ## Determining how many quarterly steps there are between 2020-10-01 and 2023-07-01.
  ## From 2023-07-01, six steps ahead will reach 2024-10-01, which is the last date with a known GDP value.
  ## Note: The first prediction date is computed by subtracting five quarters from the end of the dataset.
  ## Since the training set uses `< n_date` (excluding the observation at `n_date`), this ensures that 
  ## the final 6-step-ahead forecast lands exactly on 2024-10-01.
  start_date <- as.Date("2020-10-01")
  reference_date <- country_df$date[nrow(country_df)] %m-% months(3 * 5)  ## subtract 5 quarters 
  k <- interval(start_date, reference_date) %/% months(3) ## number of prediction rounds
  n <- k

  ## Creating empty data frames to store the predictions:
  arma_df <-  data.frame()
  trig_df <-  data.frame()
  arima_df <-  data.frame()
  
  ## Creating a data frame for the prediction number (1 to 6) - will make it easier to later filter out the 6th prediction:
  pr_nr <-  seq(1,6,1) 
  pr_nr <- pr_nr |> as.data.frame() 
  
  for (j in 1:k){
    ## Determining the first prediction date (shifts one period later each time):
    n_date <- country_df$date[nrow(country_df)-(4+n)] 
    
    ## Filtering the training set: keeping data from 1995 until (but not including) the prediction date:
    df <- country_df |> filter(date < as.Date(n_date) )
    
    
    ## Estimating the models:
    stsm_trig <- stsm_estimate(df, freq=4, seasons=F,det_drift=F, det_cycle = F, det_trend = F, decomp="trend-cycle", cycle = "trig", multiplicative = multp)
    
    stsm_arma <- stsm_estimate(df, freq=4, seasons=F,det_drift=F, det_cycle = F, det_trend = F, decomp="trend-cycle", cycle = "arma", multiplicative = multp)
    
    arima_model <- auto.arima(ts(log(df$gdp), frequency = 4, start =  c(1995, 1)), seasonal = F) 
    
  
    
    ## Finding the predictions for the next six periods:
    stsm_fitted_trig <- stsm_forecast(model=stsm_trig, y=df, n.ahead = 6) 
    stsm_fitted_arma <- stsm_forecast(model=stsm_arma, y=df, n.ahead = 6) 
    pred_arima <- forecast(arima_model, 6)
    
    
    
    ## Since stsm_forecast returns both fitted values and predictions, we need to extract only the prediction part
    trig_pred <- stsm_fitted_trig |> tail(n=6)  # Selecting the last 6 values (the actual forecasts)
    
    ## Adding the prediction starting date and the loop index to the results
    ## This helps verify the correctness of the loop and makes it easier to visually distinguish between forecasts
    trig_pred <- trig_pred |> mutate(pred_date = as.Date(n_date), k = j)
 
    
    arma_pred <- stsm_fitted_arma |> tail(n=6)
    arma_pred <- arma_pred |> mutate(pred_date = as.Date(n_date), k = j)
    
    
    ## Extracting the prediction dates (they are the same across all models):
    pred_dates <- trig_pred$date 
    ## These dates are added to the ARIMA forecasts, which do not automatically include a date column:
    pred_arima <- as.data.frame(pred_arima) |> mutate(pred_date = as.Date(n_date), k = j) 
    pred_arima <- cbind(pred_arima, date=as.Date(pred_dates))
    
    
    
    ## Combining results into a single data frame
    trig_df <- rbind(trig_df, cbind(trig_pred, pr_nr))
    arma_df <- rbind(arma_df, cbind(arma_pred, pr_nr)) 
    arima_df <- rbind(arima_df, cbind(pred_arima, pr_nr))
    
    ## Reduce n by one
    n = n-1
  }
  
  return(list(trig_df=trig_df,arma_df=arma_df,arima_df=arima_df ))
}
```



The purpose of this function is to calculate predictions for all countries and compile them into a single data frame.
```{r}
fun_all_res <- function(countries, multp){
  ## Empty data frames to store all results:
  trig_all_count <- data.frame()
  arma_all_count <- data.frame()
  arima_all_count <- data.frame()
  
  for (j in 1:length(countries)){
    ## Computing the results:
    result <- prediction_function(get(countries[j]), multp)
    
    ## Separating the results by method:
    trig_df <- result$trig_df
    arma_df <- result$arma_df
    arima_df <- result$arima_df
    
    ## Adding the corresponding country:
    trig_df <- trig_df |> mutate(country=countries[j])
    arma_df <- arma_df |> mutate(country=countries[j])
    arima_df <- arima_df |> mutate(country=countries[j])
  
    ## Combining results into a single data frame
    trig_all_count <- rbind(trig_all_count,trig_df)
    arma_all_count <- rbind(arma_all_count,arma_df)
    arima_all_count <- rbind(arima_all_count,arima_df)
    
    
  }
  return(list(trig_all_count=trig_all_count, arma_all_count=arma_all_count, arima_all_count=arima_all_count))
}

```



The purpose of this function is to calculate RMSE and MAPE in two ways: 1. using only the sixth prediction, and 2. using all six predictions.
```{r}
fun_rmse_mape <- function(results_all, countries){
  rmse_all <- data.frame() 
  mape_all <- data.frame() 
  for (j in 1:length(countries)){
    ## Separating the results based on methods and filtering out a specific country
    result_trig <- results_all$trig_all_count |> filter(country == countries[j])
    result_arma <- results_all$arma_all_count |> filter(country == countries[j])
    result_arima <- results_all$arima_all_count |> filter(country == countries[j])
    
    ## Selecting the 6th prediction:
    trig6 <- result_trig |> filter(pr_nr == 6)
    arma6 <- result_arma |> filter(pr_nr == 6)
    arima6 <- result_arima |> filter(pr_nr == 6) |> dplyr::select(pred_date, pred_arima =`Point Forecast`)
    exact_gdp <-  get(countries[j]) |> tail(n=nrow(arma6)) ## Selecting the exact gdp values
    

    ## RMSE for the 6th prediction step:
    rmse_trig_6 <- sqrt(mean((exact_gdp$gdp-trig6$fitted)**2))
    rmse_arma_6 <- sqrt(mean((exact_gdp$gdp-arma6$fitted)**2))
    rmse_arima_6 <- sqrt(mean((exact_gdp$gdp-exp(arima6$pred_arima) )**2))
    
    
    ## MAPE for the 6th prediction step:
    mape_trig_6 <- (mean(abs(exact_gdp$gdp-trig6$fitted)/exact_gdp$gdp))
    mape_arma_6 <-  (mean(abs(exact_gdp$gdp-arma6$fitted)/exact_gdp$gdp))
    mape_arima_6 <- (mean(abs(exact_gdp$gdp-exp(arima6$pred_arima) )/exact_gdp$gdp))
        
    
    ## ---
    
    ## Calculating how many prediction dates there are to extract the exact same number of observed GDP values
    ## This ensures that we only select the values needed for evaluating accuracy
    k <- result_trig |> distinct(result_trig |> dplyr::select(date)) |> count()
    exact_gdp <-  get(countries[j]) |> tail(n=k[[1]]) 
    
    trig <- merge(result_trig, exact_gdp, by = "date")
    arma <- merge(result_arma, exact_gdp, by = "date") 
    arima <- merge(result_arima, exact_gdp, by = "date") 
    
    
    ## RMSE (across all 6 steps):
    rmse_trig <- trig |> summarise(rmse = sqrt(mean((gdp-fitted)**2))) 
    rmse_arma <- arma |> summarise(rmse = sqrt(mean((gdp-fitted)**2)))
    rmse_arima <- arima |> summarise(rmse = sqrt(mean((gdp-exp(`Point Forecast`) )**2))) 
  
    
    ## MAPE (across all 6 steps):
    mape_trig <- trig |> summarise(mape = (mean(abs(gdp-fitted)/gdp)))
    mape_arma <- arma |> summarise(mape = (mean(abs(gdp-fitted)/gdp)))
    mape_arima <- arima |> summarise(mape = (mean(abs(gdp-exp(`Point Forecast`) )/gdp)))
     
    
    ## Combining the results into one single row that will be included in the overall results table:
    row_rmse <- c(countries[j],rmse_trig_6, rmse_arma_6, rmse_arima_6,  rmse_trig[[1]], rmse_arma[[1]], rmse_arima[[1]])
    row_mape <- c(countries[j], mape_trig_6, mape_arma_6, mape_arima_6, mape_trig[[1]], mape_arma[[1]], mape_arima[[1]]) 
    
    ## Combining results into a single data frame
    rmse_all <- rbind(rmse_all, row_rmse)
    mape_all <- rbind(mape_all, row_mape)
    
  }
  colnames(rmse_all) <- c("country", "rmse_trig_6", "rmse_arma_6", "rmse_arima_6", 
                           "rmse_trig" ,  "rmse_arma",  "rmse_arima")
  
  colnames(mape_all) <- c("country", "mape_trig_6", "mape_arma_6", "mape_arima_6",
                          "mape_trig" ,  "mape_arma",  "mape_arima")
  return(list(rmse_all=rmse_all, mape_all=mape_all))
}

```



## RUN THIS
```{r}
results_all <- fun_all_res(c("estonia", "latvia", "lithuania", "germany", "sweden", "finland"), TRUE)
```


```{r}
rmse_mape_results<- fun_rmse_mape(results_all,  c("estonia", "latvia", "lithuania", "germany", "sweden", "finland")) 

rmse_mape_results
```


```{r}
#writexl::write_xlsx(rmse_mape_results, "mltp_rmse_mape_results.xlsx")
```



***


***


#### Plots Used in the Thesis:

***
## Economic growth plot:

##### Calculating the economic growth:
```{r}
## NB! Once this code has been run
estonia <- estonia |> mutate(growth=100*(gdp - lag(gdp, 4)) / lag(gdp, 4))
finland <- finland |> mutate(growth=100*(gdp - lag(gdp, 4)) / lag(gdp, 4))
latvia <- latvia  |> mutate(growth=100*(gdp - lag(gdp, 4)) / lag(gdp, 4))
sweden <- sweden  |> mutate(growth=100*(gdp - lag(gdp, 4)) / lag(gdp, 4))
germany <- germany |> mutate(growth=100*(gdp - lag(gdp, 4)) / lag(gdp, 4))
lithuania <- lithuania  |> mutate(growth=100*(gdp - lag(gdp, 4)) / lag(gdp, 4))
```


#### Plots of the economic growth:
```{r message=FALSE, warning=FALSE}
economic_growth1 <- ggplot() +
  geom_line(data = estonia,  mapping = aes(date, growth, color = "Estonia")) +
  geom_line(data = germany,  mapping = aes(date, growth, color = "Germany")) +
  geom_line(data = sweden,   mapping = aes(date, growth, color = "Sweden")) +
  geom_line(data = finland,  mapping = aes(date, growth, color = "Finland")) +
  geom_hline(yintercept = 0, color = "darkgray") + 
  scale_color_manual(values = c("Estonia" = "#0072B2", "Germany" = "black", "Sweden" = "#F5C800", "Finland" = "#33B7D8")) +
  labs(x = "Date", y = "Annual GDP growth (%)", color = "Country") +
  theme_minimal()  

ggsave("economic_growth1.jpg", economic_growth1, width=7.5, height=5)


economic_growth2 <- ggplot() +
  geom_line(data = estonia,  mapping = aes(date, growth, color = "Estonia")) +
  geom_line(data = latvia,  mapping = aes(date, growth, color = "Latvia")) +
  geom_line(data = lithuania,   mapping = aes(date, growth, color = "Lithuania")) +
  scale_color_manual(values =  c("Estonia" = "#0072B2", "Latvia" = "#A40000", "Lithuania" = "#E69F00")) +
  geom_hline(yintercept = 0, color = "darkgray") + 
  labs(x = "Date", y = "Annual GDP growth (%)", color = "Country") +
  theme_minimal()  

ggsave("economic_growth2.jpg", economic_growth2, width=7.5, height=5)

```


***




#### Comparing methods:
```{r}
plot_fun_methods <- function(df, CN){ 
  df_hp <- hpfilter(log(df$gdp),freq=1600,type=c("lambda"), drift = T)
  df_trig_fitted <- trig_models |> filter(country==CN)
  df_arma_fitted <- arma_models |> filter(country==CN)

  data_plot <- data.frame(
    date = df$date,
    hp_cycle = df_hp$cycle,
    trig_cycle = (df_trig_fitted$cycle-1), #/df_trig_fitted$observed
    arma_cycle = (df_arma_fitted$cycle-1)) #/df_arma_fitted$observed
  colnames(data_plot) <- c("date", "HP cycle", "trig cycle", "ARMA cycle")
  
  data_long <- data_plot |> 
    tidyr::pivot_longer(cols = -date, names_to = "method", values_to = "cycle")

  name <- paste0(toupper(substring(CN, 1, 1)), substring(CN, 2))
  plot <- ggplot(data_long, aes(x = date, y = cycle, color = method)) +
   geom_line() +
   geom_hline(yintercept = 0, color = "darkgray") + 
   scale_color_manual(values = c("HP cycle" = "#1F90B7", "trig cycle" = "#11B90B", "ARMA cycle" = "#e9091d")) +
   labs(x = "Date", y = "Cycle", color = "Method:") + ggtitle(name)+
   theme_minimal()+
  theme(legend.position = "bottom")
  
  return(plot)
  
}


plot_ee <- plot_fun_methods(estonia,"estonia") 
ggsave("mltp_estonia_cycle_methods.jpg", plot_ee, width=7.5, height=5)

plot_de <- plot_fun_methods(germany,"germany") 
ggsave("mltp_germany_cycle_methods.jpg", plot_de, width=7.5, height=5)

plot_se <- plot_fun_methods(sweden,"sweden") 
ggsave("mltp_sweden_cycle_methods.jpg", plot_se, width=7.5, height=5)

plot_lv <- plot_fun_methods(latvia,"latvia") 
ggsave("mltp_latvia_cycle_methods.jpg", plot_lv, width=7.5, height=5)

plot_lt <- plot_fun_methods(lithuania,"lithuania") 
ggsave("mltp_lithuania_cycle_methods.jpg", plot_lt, width=7.5, height=5)

plot_fi <- plot_fun_methods(finland,"finland") 
ggsave("mltp_finland_cycle_methods.jpg", plot_fi, width=7.5, height=5)
```




#### Comparing results of different countries (trig):
```{r}
trig_ee <- trig_models |> filter(country=="estonia")
trig_lv <- trig_models |> filter(country=="latvia")
trig_lt <- trig_models |> filter(country=="lithuania")
trig_de <- trig_models |> filter(country=="germany")
trig_fi <- trig_models |> filter(country=="finland")
trig_se <- trig_models |> filter(country=="sweden")


trig_cycle_baltics <- ggplot() +
  geom_line(data = trig_ee,  mapping = aes(date, cycle-1, color = "Estonia")) +
  geom_line(data = trig_lv,  mapping = aes(date, cycle-1, color = "Latvia")) +
  geom_line(data = trig_lt,   mapping = aes(date, cycle-1, color = "Lithuania")) +
  geom_hline(yintercept = 0, color = "darkgray") +
  scale_color_manual(values = c("Estonia" = "#0072B2", "Latvia" = "#A40000", "Lithuania" = "#E69F00")) +
  labs(x = "Date", y = "Cycle", color = "Country:") +
  theme_minimal() +
  theme(legend.position = "bottom") 

ggsave("mltp_trig_cycle_baltics.jpg", trig_cycle_baltics, width=7.5, height=5)

trig_cycle_others <- ggplot() +
  geom_line(data = trig_ee,  mapping = aes(date, cycle-1, color = "Estonia")) +
  geom_line(data = trig_de,  mapping = aes(date, cycle-1, color = "Germany")) +
  geom_line(data = trig_se,   mapping = aes(date, cycle-1, color = "Sweden"))+
  geom_line(data = trig_fi,   mapping = aes(date, cycle-1, color = "Finland"))+
  geom_hline(yintercept = 0, color = "darkgray") +
  scale_color_manual(values = c("Estonia" = "#0072B2", "Germany" = "black", "Sweden" = "#F5C800", "Finland" = "#33B7D8")) + 
  labs(x = "Date", y = "Cycle ", color = "Country:") +
  theme_minimal() +
  theme(legend.position = "bottom") 

ggsave("mltp_trig_cycle_others.jpg", trig_cycle_others, width=7.5, height=5)

#------ Trend -------#



trig_trend_baltics <- ggplot() +
  geom_line(data = trig_ee,  mapping = aes(date, trend/observed[1], color = "Estonia")) +
  geom_line(data = trig_lv,  mapping = aes(date, trend/observed[1], color = "Latvia")) +
  geom_line(data = trig_lt,   mapping = aes(date, trend/observed[1], color = "Lithuania")) +
  scale_color_manual(values = c("Estonia" = "#0072B2", "Latvia" = "#A40000", "Lithuania" = "#E69F00")) +
  labs(x = "Date", y = "Trend / first observed value", color = "Country:") +
  theme_minimal() +
  theme(legend.position = "bottom") 


ggsave("mltp_trig_trend_baltics.jpg", trig_trend_baltics, width=7.5, height=5)

trig_trend_others <-ggplot() +
  geom_line(data = trig_ee,  mapping = aes(date, trend/observed[1], color = "Estonia")) +
  geom_line(data = trig_de,  mapping = aes(date, trend/observed[1], color = "Germany")) +
  geom_line(data = trig_se,   mapping = aes(date, trend/observed[1], color = "Sweden"))+
  geom_line(data = trig_fi,   mapping = aes(date, trend/observed[1], color = "Finland"))+
  scale_color_manual(values = c("Estonia" = "#0072B2", "Germany" = "black", "Sweden" = "#F5C800", "Finland" = "#33B7D8")) +
  labs(x = "Date", y = "Trend / first observed value", color = "Country:") +
  theme_minimal()  +
  theme(legend.position = "bottom")


ggsave("mltp_trig_trend_others.jpg", trig_trend_others, width=7.5, height=5)

```


#### Comparing results of different countries (arma):
```{r}
arma_ee <- arma_models |> filter(country=="estonia")
arma_lv <- arma_models |> filter(country=="latvia")
arma_lt <- arma_models |> filter(country=="lithuania")
arma_de <- arma_models |> filter(country=="germany")
arma_fi <- arma_models |> filter(country=="finland")
arma_se <- arma_models |> filter(country=="sweden")


arma_cycle_baltics <- ggplot() +
  geom_line(data = arma_ee,  mapping = aes(date, cycle-1, color = "Estonia")) +
  geom_line(data = arma_lv,  mapping = aes(date, cycle-1, color = "Latvia")) +
  geom_line(data = arma_lt,   mapping = aes(date, cycle-1, color = "Lithuania")) +
  geom_hline(yintercept = 0, color = "darkgray") +
  scale_color_manual(values = c("Estonia" = "#0072B2", "Latvia" = "#A40000", "Lithuania" = "#E69F00")) +
  labs(x = "Date", y = "Cycle", color = "Country:") +
  theme_minimal()  +
  theme(legend.position = "bottom")

ggsave("mltp_arma_cycle_baltics.jpg", arma_cycle_baltics, width=7.5, height=5)

arma_cycle_others <- ggplot() +
  geom_line(data = arma_ee,  mapping = aes(date, cycle-1, color = "Estonia")) +
  geom_line(data = arma_de,  mapping = aes(date, cycle-1, color = "Germany")) +
  geom_line(data = arma_se,   mapping = aes(date, cycle-1, color = "Sweden"))+
  geom_line(data = arma_fi,   mapping = aes(date, cycle-1, color = "Finland"))+
  geom_hline(yintercept = 0, color = "darkgray") +
  scale_color_manual(values = c("Estonia" = "#0072B2", "Germany" = "black", "Sweden" = "#F5C800", "Finland" = "#33B7D8")) +
  labs(x = "Date", y = "Cycle", color = "Country:") +
  theme_minimal()  +
  theme(legend.position = "bottom")

ggsave("mltp_arma_cycle_others.jpg", arma_cycle_others, width=7.5, height=5)


#------ Trend -------#


arma_trend_baltics <- ggplot() +
  geom_line(data = arma_ee,  mapping = aes(date, trend/observed[1], color = "Estonia")) +
  geom_line(data = arma_lv,  mapping = aes(date, trend/observed[1], color = "Latvia")) +
  geom_line(data = arma_lt,   mapping = aes(date, trend/observed[1], color = "Lithuania")) +
  scale_color_manual(values = c("Estonia" = "#0072B2", "Latvia" = "#A40000", "Lithuania" = "#E69F00")) +
  labs(x = "Date", y = "Trend / first observed value", color = "Country:") +
  theme_minimal() +
  theme(legend.position = "bottom") 

ggsave("mltp_arma_trend_baltics.jpg", arma_trend_baltics, width=7.5, height=5)

arma_trend_others <- ggplot() +
  geom_line(data = arma_ee,  mapping = aes(date, trend/observed[1], color = "Estonia")) +
  geom_line(data = arma_de,  mapping = aes(date, trend/observed[1], color = "Germany")) +
  geom_line(data = arma_se,   mapping = aes(date, trend/observed[1], color = "Sweden"))+
  geom_line(data = arma_fi,   mapping = aes(date, trend/observed[1], color = "Finland"))+
  scale_color_manual(values = c("Estonia" = "#0072B2", "Germany" = "black", "Sweden" = "#F5C800", "Finland" = "#33B7D8")) +
  labs(x = "Date", y = "Trend / first observed value", color = "Country:") +
  theme_minimal()  +
  theme(legend.position = "bottom")  

ggsave("mltp_arma_trend_others.jpg", arma_trend_others, width=7.5, height=5)

```



#### Forecasted vs observed:
```{r}
plot_forecast <- function(CN){ 
  result_trig <- results_all$trig_all_count |> filter(country == CN)
  result_arma <- results_all$arma_all_count |> filter(country == CN)
  result_arima <- results_all$arima_all_count |> filter(country == CN)
    
  trig6 <- result_trig |> filter(pr_nr == 6)
  arma6 <- result_arma |> filter(pr_nr == 6)
  arima6 <- result_arima |> filter(pr_nr == 6) 
  exact_gdp <-  get(CN) |> tail(n=nrow(arma6))
  
  name <- paste0(toupper(substring(CN, 1, 1)), substring(CN, 2))
   plot <- ggplot() +
  geom_line(data = exact_gdp,  mapping = aes(date, gdp, color = "GDP")) +
  geom_point(data = exact_gdp,  mapping = aes(date, gdp, color = "GDP")) + 
  geom_line(data = trig6,  mapping = aes(date, fitted , color = "Trig prediction")) +
  geom_point(data = trig6,  mapping = aes(date, fitted , color = "Trig prediction"))+
  geom_line(data = arma6,   mapping = aes(date, fitted , color = "ARMA prediction")) +
  geom_point(data = arma6,   mapping = aes(date, fitted , color = "ARMA prediction"))+
  geom_line(data = arima6,   mapping = aes(date, exp(`Point Forecast`) , color = "ARIMA prediction")) +
  geom_point(data = arima6,   mapping = aes(date, exp(`Point Forecast`) , color = "ARIMA prediction"))+
  scale_color_manual(values = c("GDP" = "#000000", "Trig prediction" = "#11B90B", "ARMA prediction" = "#e9091d", "ARIMA prediction" = "#1F90B7"), breaks = c("GDP", "Trig prediction", "ARMA prediction", "ARIMA prediction")) +
  labs(x = "Date", y = "GDP", color = "Forecast:") + ggtitle(name)+
  theme_minimal()   +
  theme(legend.position = "bottom")
  
  return(plot)
}  

pred_ee <- plot_forecast("estonia")
ggsave("mltp_forecast_estonia.jpg", pred_ee, width=7.5, height=5)

pred_de <- plot_forecast("germany")
ggsave("mltp_forecast_germany.jpg", pred_de, width=7.5, height=5)

pred_lv <- plot_forecast("latvia")
ggsave("mltp_forecast_latvia.jpg", pred_lv, width=7.5, height=5)

pred_lt <- plot_forecast("lithuania")
ggsave("mltp_forecast_lithuania.jpg", pred_lt, width=7.5, height=5)

pred_se <- plot_forecast("sweden")
ggsave("mltp_forecast_sweden.jpg", pred_se, width=7.5, height=5)

pred_fi <- plot_forecast("finland")
ggsave("mltp_forecast_finland.jpg", pred_fi, width=7.5, height=5)
```




***

### Model theoretical form:
```{r}
CN <- "estonia"
coef_trig <- trig_model_forms |> dplyr::filter(country==CN) |> dplyr::select(coef)
coef_trig[[1]]  

CN <- "estonia"
coef_arma <- arma_model_forms |> dplyr::filter(country==CN) |> dplyr::select(coef)
coef_arma[[1]]  
```



Y_t = T_t + C_t + e_t, e_t ~ NID(0, sig_e^2)      

T_t = T_{t-1} + D_{t-1} + e_t, e_t ~ NID(0, sig_t^2)

D_t = d + phi_d*D_{t-1} + n_t;  n_t ~ NID(0, sig_d^2)

C_t = ... + k_t; k_t ~ NID(0, sig_c^2)



Y_t = T_t + C_t + e_t               -- Main equation
T_t = T_{t-1} + D_{t-1} + e_t       -- Level equation
D_t = d + phi_d*D_{t-1} + n_t       -- Slope equation

C_t = ... + k_t                     -- Cycle equation (trig/ARMA formulation)

***

***