gdp_w_gold.Rmd

---
title: "GDP y salario en Oro"
output: html_notebook
---

Me quedo con el __GDP__ y el __Salario__, lo normalizo para que este expresado en __oro__

```{r message=FALSE, warning=FALSE}
library(tidyverse)
library(lubridate)
library(forecast)
library(magrittr)
library(broom)
library(sweep)
library(ggthemes)

```


```{r load and basic changes}
gold <- read_csv("data/GOLD_1791-2018.csv",skip = 3)
gdp <- read_csv("data/USGDP_1790-2018.csv", skip=2)
wage <- read_csv("data/USWAGE_1774-2018.csv", skip=3)

gdp_in_gold <- left_join(gold, gdp, by = "Year") %>% 
  transmute(value = `Nominal GDP per capita (current dollars)`/`New York Market Price (U.S. dollars per fine ounce)`,
         year = parse_date_time(Year,"y")) 
wg_in_gold <- wage %>% 
  #filter(Year>=1900) %>% 
  left_join(gold, gdp, by = "Year") %>% 
  transmute(value =`Production Workers Hourly Compensation (nominal dollars)`/`New York Market Price (U.S. dollars per fine ounce)`,
            year = parse_date_time(Year,"y")) %>%
  na.omit()


#lo paso a serie de tiempo
wg_in_gold$value <- ts(wg_in_gold$value,start = min(year(wg_in_gold$year)))
gdp_in_gold$value <- ts(gdp_in_gold$value,start = min(year(gdp_in_gold$year)))

```


```{r}
range(gdp_in_gold$year)
range(wg_in_gold$year)
```

Agrupo las series antes y después de 1971

```{r}

df <- wg_in_gold %>% 
  mutate(tipo = "wage") %>% 
  bind_rows(gdp_in_gold %>% 
              mutate(tipo = "gdp")) %>% 
    mutate(group=case_when(year< parse_date_time(1971,orders = 'y')~"Antes del 71'",
                         TRUE ~ "Despues del 71'"))
```


## Autocorrelación, todo el período

```{r}
ggAcf(diff(df %>% 
             filter(tipo=="gdp") %$% value ))+
  labs(title= "GDP. Todo el período")+
  theme_minimal()

ggsave("plots/Autocorrelacion_gdp.png",dpi=300)

ggAcf(diff(df %>% 
             filter(tipo=="wage") %$% value ))+
  labs(title= "Salario. Todo el período")+
  theme_minimal()

ggsave("plots/Autocorrelacion_wg.png",dpi=300)

```


#### GDP
la serie en diferencias tiene una autocorrelación importante con el primer lag (tipo random walk) y una correlación negativa con 10 años atrás.


#### Salarios
Con la serie entera, además del período anterior, tiene una relación positiva con lo que sucede 6 años atrás (ídem GDP)



## Antes y después del 71


```{r}
diferencio <- function(df){diff(df$value)}

nested_df <- df %>% 
  group_by(tipo,group) %>%
  nest() %>% 
  mutate(diferencias  = map(data,diferencio))
  
nested_df
```


```{r}
nested_df %>% 
  filter(tipo == "gdp", group=="Antes del 71'") %$%
  ggAcf(diferencias[[1]])+
  labs(title = "GDP. Antes de 1971")+  
  theme_minimal()


ggsave("plots/Autocorrelacion_gdp_pre71.png",dpi=300)

nested_df %>% 
  filter(tipo == "gdp", group=="Despues del 71'") %$%
  ggAcf(diferencias[[1]])+
  labs(title = "GDP. Después de 1971")+
  theme_minimal()

ggsave("plots/Autocorrelacion_gdp_post71.png",dpi=300)

nested_df %>% 
  filter(tipo == "wage", group=="Antes del 71'") %$%
  ggAcf(diferencias[[1]])+
  labs(title = "Salario. Antes de 1971")+
  theme_minimal()


ggsave("plots/Autocorrelacion_wg_pre71.png",dpi=300)

nested_df %>% 
  filter(tipo == "wage", group=="Despues del 71'") %$%
  ggAcf(diferencias[[1]])+
  labs(title = "Salario. Después de 1971")+
  theme_minimal()


ggsave("plots/Autocorrelacion_wg_post71.png",dpi=300)

```

__GDP__: Es como si con el patrón oro se pudiera ver una autocorrealción para el ciclo corto (3-5 años, _real business cycle_) pero a partir de los 70 es un Random Walk

__Salario__:
- Antes de los 70' no tiene autocorrelación significativa con ningún período (¿aleatorio puro??)
- Post 70' es un Random Walk


## Arima

```{r}

fit.arima <- function(df,order){
  Arima(df$value, order = order)
}


nested_df <- nested_df %>% 
  group_by(tipo, group) %>% 
  mutate(fitted_arima_2_1_0 = map(data,fit.arima, order=c(2,1,0)),
         fitted_arima_5_1_0 = map(data,fit.arima, order=c(5,1,0)),
         fitted_arima_10_1_0 = map(data,fit.arima, order=c(10,1,0)),
         fitted_arima_1_1_1 = map(data,fit.arima, order=c(1,1,1)))

nested_df <- nested_df %>% 
  mutate(tdy_2_1_0 = map(fitted_arima_2_1_0,sw_tidy),
         tdy_5_1_0 = map(fitted_arima_5_1_0,sw_tidy),
         tdy_10_1_0 = map(fitted_arima_10_1_0,sw_tidy),
         tdy_1_1_1 = map(fitted_arima_1_1_1,sw_tidy))
  
nested_df %>% 
unnest(tdy_2_1_0, .drop=TRUE)
nested_df %>% 
unnest(tdy_5_1_0, .drop=TRUE)
nested_df %>% 
unnest(tdy_10_1_0, .drop=TRUE)
nested_df %>% 
unnest(tdy_1_1_1, .drop=TRUE)

```

```{r}
nested_df %>% 
unnest(tdy_5_1_0, .drop=TRUE) %>% 
  ggplot(.,aes(term,estimate, fill= term))+
  geom_col()+
  facet_wrap(tipo~group)+
  theme_minimal()+
  theme(legend.position = "none")+
  labs(title = "AR 5, Integrado de orden 1")

```



```{r}
glances <- nested_df %>%
  mutate(g1 = map(fitted_arima_2_1_0,sw_glance), 
         g2 = map(fitted_arima_5_1_0,sw_glance),
         g3 = map(fitted_arima_10_1_0,sw_glance),
         g4 = map(fitted_arima_1_1_1,sw_glance)) %>%
  select(tipo, group, g1, g2,g3,g4) %>% 
  gather(model, glance, 3:6) %>% 
  unnest(glance, .drop = TRUE)


glances
```

El __sigma__ es la raíz de la variabilidad estimada de los residuos. Puedo comparar por serie

[MASE](https://en.wikipedia.org/wiki/Mean_absolute_scaled_error)


```{r}
ggplot(glances,aes(x=model.desc,y= sigma, color =group, size=MASE))+
geom_point()+
facet_wrap(tipo~., scales="free")+
theme(axis.text.x = element_text(angle=35, vjust = 0.5))

```

- Tanto para el gdp como para wage, la variabilidad de los residuos estimados es mayor Después de 1971 que antes. 
- En términos del poder predictivo de los modelos (MASE), el gdp antes del 71' pareciera ser lo menos predecible y luego el salario después del 71'. Para el gdp después del 71' el salario en el primer período, pareciera que el ARIMA(10,1,0) es el que logra una mejor capacidad de predicción. 

### chequeo unit roots

ARIMA(2,1,0)

```{r}

nested_df %>% 
  filter(tipo == "gdp", group=="Antes del 71'") %$%
  autoplot(fitted_arima_2_1_0[[1]])+
  labs(title = "GDP. Antes de 1971")

nested_df %>% 
  filter(tipo == "gdp", group=="Despues del 71'") %$%
  autoplot(fitted_arima_2_1_0[[1]])+
  labs(title = "GDP. Después de 1971")

nested_df %>% 
  filter(tipo == "wage", group=="Antes del 71'") %$%
  autoplot(fitted_arima_2_1_0[[1]])+
  labs(title = "Salario. Antes de 1971")

nested_df %>% 
  filter(tipo == "wage", group=="Despues del 71'") %$%
  autoplot(fitted_arima_2_1_0[[1]])+
  labs(title = "Salario. Después de 1971")

```



ARIMA(5,1,0)

```{r}

nested_df %>% 
  filter(tipo == "gdp", group=="Antes del 71'") %$%
  autoplot(fitted_arima_5_1_0[[1]])+
  labs(title = "GDP. Antes de 1971")

nested_df %>% 
  filter(tipo == "gdp", group=="Despues del 71'") %$%
  autoplot(fitted_arima_5_1_0[[1]])+
  labs(title = "GDP. Después de 1971")

nested_df %>% 
  filter(tipo == "wage", group=="Antes del 71'") %$%
  autoplot(fitted_arima_5_1_0[[1]])+
  labs(title = "Salario. Antes de 1971")

nested_df %>% 
  filter(tipo == "wage", group=="Despues del 71'") %$%
  autoplot(fitted_arima_5_1_0[[1]])+
  labs(title = "Salario. Después de 1971")

```


ARIMA(10,1,0)

```{r}

nested_df %>% 
  filter(tipo == "gdp", group=="Antes del 71'") %$%
  autoplot(fitted_arima_10_1_0[[1]])+
  labs(title = "GDP. Antes de 1971")

nested_df %>% 
  filter(tipo == "gdp", group=="Despues del 71'") %$%
  autoplot(fitted_arima_10_1_0[[1]])+
  labs(title = "GDP. Después de 1971")

nested_df %>% 
  filter(tipo == "wage", group=="Antes del 71'") %$%
  autoplot(fitted_arima_10_1_0[[1]])+
  labs(title = "Salario. Antes de 1971")

nested_df %>% 
  filter(tipo == "wage", group=="Despues del 71'") %$%
  autoplot(fitted_arima_10_1_0[[1]])+
  labs(title = "Salario. Después de 1971")

```


ARIMA(1,1,1)

```{r}

nested_df %>% 
  filter(tipo == "gdp", group=="Antes del 71'") %$%
  autoplot(fitted_arima_1_1_1[[1]])+
  labs(title = "GDP. Antes de 1971")

nested_df %>% 
  filter(tipo == "gdp", group=="Despues del 71'") %$%
  autoplot(fitted_arima_1_1_1[[1]])+
  labs(title = "GDP. Después de 1971")

nested_df %>% 
  filter(tipo == "wage", group=="Antes del 71'") %$%
  autoplot(fitted_arima_1_1_1[[1]])+
  labs(title = "Salario. Antes de 1971")

nested_df %>% 
  filter(tipo == "wage", group=="Despues del 71'") %$%
  autoplot(fitted_arima_1_1_1[[1]])+
  labs(title = "Salario. Después de 1971")

```


TODO adentro del circulo unitario


quiero ver cómo quedan los suavizados


```{r}

nested_df %>%
  mutate(augment = map(fitted_arima_2_1_0,sw_augment)) %>% 
  unnest(augment, .drop = TRUE)

nested_df %>%
  mutate(augment = map(fitted_arima_5_1_0,sw_augment)) %>% 
  unnest(augment, .drop = TRUE)

nested_df %>%
  mutate(augment = map(fitted_arima_10_1_0,sw_augment)) %>% 
  unnest(augment, .drop = TRUE)

nested_df %>%
  mutate(augment = map(fitted_arima_1_1_1,sw_augment)) %>% 
  unnest(augment, .drop = TRUE)

```

```{r fig.height=10,fig.width=10}
nested_df %>%
  mutate(augment_2_1_0 = map(fitted_arima_2_1_0,sw_augment),
         augment_5_1_0 = map(fitted_arima_5_1_0,sw_augment),
         augment_10_1_0 = map(fitted_arima_10_1_0,sw_augment),
         augment_1_1_1 = map(fitted_arima_1_1_1,sw_augment)) %>% 
  select(tipo, group, data, augment_2_1_0, augment_5_1_0,  augment_10_1_0, augment_1_1_1) %>%
   gather(modelo, augment, 4:ncol(.)) %>% 
  unnest() %>% 
  ggplot(.,aes(year,.fitted, color = modelo))+
  geom_line()+
  facet_wrap(tipo~group, scales = "free")+
  theme_tufte()+
  theme(legend.position = "bottom")
```


¯\\_(ツ)_/¯


## Auto.Arima


```{r}

nested_df <- df %>% 
  group_by(tipo,group) %>%
  nest() 

fit_auto_arima <- function(df){
  auto.arima(df$value)
}


nested_df <- nested_df %>% 
  group_by(tipo, group) %>% 
  mutate(fitted_auto_arima = map(data,fit_auto_arima),
         tidy = map(fitted_auto_arima,sw_tidy),
         glance = map(fitted_auto_arima,sw_glance),
         augment = map(fitted_auto_arima,sw_augment))
```


```{r}
nested_df %>% 
unnest(glance, .drop=TRUE) %>% 
  select(tipo, group, model.desc, sigma,RMSE) %>% 
  xtable::xtable(.)


```


Antes del 71' para ambas variables el mejor modelo es un proceso integrado de orden dos, mientras que para las series después del 71' el mejor modelo no es integrado, pero tiene una media distinta a cero. 

Además, antes de 71' los modelos tienen menor desvío y RMSE


```{r}
nested_df %>% 
  filter(tipo == "gdp", group=="Antes del 71'") %$%
  autoplot(fitted_auto_arima[[1]])+
  labs(title = "GDP. Antes de 1971")

nested_df %>% 
  filter(tipo == "gdp", group=="Despues del 71'") %$%
  autoplot(fitted_auto_arima[[1]])+
  labs(title = "GDP. Después de 1971")

nested_df %>% 
  filter(tipo == "wage", group=="Antes del 71'") %$%
  autoplot(fitted_auto_arima[[1]])+
  labs(title = "Salario. Antes de 1971")

nested_df %>% 
  filter(tipo == "wage", group=="Despues del 71'") %$%
  autoplot(fitted_auto_arima[[1]])+
  labs(title = "Salario. Después de 1971")

```



```{r}
nested_df %>% 
  select(tipo, group, data, augment) %>%
  unnest() %>% 
  ggplot(.,aes(year))+
  geom_line(aes(y=.fitted,color="fitted"))+
  geom_line(aes(y=value,color="real"))+
  facet_wrap(tipo~group, scales = "free")+
  theme_tufte()+
  theme(legend.position = "bottom")

ggsave("plots/arima_groups.png", dpi = 300)
```


El modelo funciona mejor para antes de 71' que para después.


Sin quebrar las series:



```{r}

nested_df <- df %>%
  select(-group) %>% 
  group_by(tipo) %>%
  nest() 

nested_df <- nested_df %>% 
  group_by(tipo) %>% 
  mutate(fitted_auto_arima = map(data,fit_auto_arima),
         tidy = map(fitted_auto_arima,sw_tidy),
         glance = map(fitted_auto_arima,sw_glance),
         augment = map(fitted_auto_arima,sw_augment))

nested_df %>% 
  select(tipo, data, augment) %>%
  unnest() %>% 
  ggplot(.,aes(year))+
  geom_line(aes(y=.fitted,color="fitted"))+
  geom_line(aes(y=value,color="real"))+
  facet_grid(tipo~., scales = "free")+
  theme_tufte()+
  theme(legend.position = "bottom")

ggsave("plots/arima_tot.png", dpi = 300)

```


```{r}
nested_df %>% 
  unnest(glance)
```

El RMSE queda en el medio entre ambos. No pareciera mejorar la predicción general separando.