wavelets.Rmd

---
title: "Wavelets"
output: html_notebook
---


```{r setup}
library(WaveletComp)
library(tidyverse)
library(lubridate)
library(glue)
library(ggthemes)
library(cowplot)
library(scales) # to access breaks/formatting functions

```


## datos
```{r message=FALSE, warning=FALSE}
gold <- read_csv("data/GOLD_1791-2018.csv",skip = 2)
gdp <- read_csv("data/USGDP_1790-2018.csv", skip=2)
wage <- read_csv("data/USWAGE_1774-2018.csv")
ir <- read_csv("data/INTERESTRATE_1857-2018.csv", skip=1)
ipc <- read_csv("data/USCPI_1774-2018.csv", skip=3)

uk_gdp <- read_csv("data/gdp_uk_gold.csv")
```

[WaveletComp](http://www.hs-stat.com/projects/WaveletComp/WaveletComp_guided_tour.pdf)

# Modelo teórico

construyo el _modelo teórico_ de una economía cíclica con que se construye a partrir de:


![periodo y amplitud](https://www.mathsisfun.com/algebra/images/period-amplitude.svg)

- un cambio de nivel (cte)
- Una tendencia (x0)
- Un ciclo corto (de periodo de 3 años y amplitud 10)
- Un ciclo medio (de periodo de 10 años y amplitud 20)
- Un ciclo largo (de periodo de 50 años y amplitud 30)
- Ruido normal

```{r}
nn = 1000
#cambio_nivel
impulso= c(rep(50,(nn/2-1)),100,rep(50,nn/2))
x0 = ts(c(1:nn)/2)
x3 = ts(10*sin((2*pi/3)*c(1:nn)))
x10 = ts(20*sin((2*pi/10)*c(1:nn)))
x50 = ts(30*sin((2*pi/50)*c(1:nn)))
ruido <- rnorm(nn)
x = impulso+ x0+  x3 + x10 + x50+ruido

```

```{r}
df <- data_frame(periodo=1:nn,impulso=impulso,tendencia=x0, ciclo_3=x3,ciclo_10=x10, ciclo_50=x50,ruido, composicion_series=x)
df %>%
  filter(periodo %in% c(450:550)) %>% 
  gather(componente,valor,2:8,factor_key=T) %>% 
  ggplot(.,aes(periodo,valor))+
  geom_line()+
  labs(title= "Elementos de la serie teórica")+
  facet_wrap(~componente, scales="free")+
  theme_minimal()

ggsave("plots/serie_teorica.PNG")

```


```{r}
df <- data_frame(period=1:nn,impulse=impulso,trend=x0, short=x3,middle=x10, long=x50,noise=ruido, composite_series=x)

df %>%
  filter(period %in% c(450:550)) %>% 
  gather(componente,valor,2:8,factor_key=T) %>% 
  ggplot(.,aes(period,valor))+
  geom_line()+
  labs(y='value')+
  facet_wrap(~componente, scales="free")+
  theme_minimal()

ggsave("plots/serie_teorica_en.PNG")

```







### Teoria ciclo


```{r}
df <- data_frame(periodo=1:nn,impulso=impulso,tendencia=x0, ciclo_3=x3,ciclo_10=x10, ciclo_50=x50,ruido, composicion_series=x)

df %>% 
  filter(periodo %in% 0:100) %>% 
ggplot(.,aes(x=periodo, as.numeric(ciclo_50)))+
  geom_line()+
  geom_segment(aes(x=25, xend=75, y=0, yend=0,color = "periodo"))+
  geom_text(aes(x=60, y=2, label="periodo", color="periodo"), size=6)+
  geom_segment(aes(x=13, xend=13, y=-30, yend=30,color = "amplitud"))+
  geom_text(aes(x=10, y=2, label="amplitud", color="amplitud"),angle=90, size=6)+
  scale_color_gdocs()+
  theme_void()+
  theme(legend.position = "none")

ggsave("plots/ciclo.png")
```



```{r}
df <- data_frame(period=1:nn,impulse=impulso,trend=x0, short=x3,middle=x10, long=x50,noise=ruido, composite_series=x)

df %>% 
  filter(period %in% 0:100) %>% 
ggplot(.,aes(x=period, as.numeric(long)))+
  geom_line()+
  geom_segment(aes(x=25, xend=75, y=0, yend=0,color = "period"))+
  geom_text(aes(x=60, y=2, label="period", color="period"), size=6)+
  geom_segment(aes(x=13, xend=13, y=-30, yend=30,color = "amplitude"))+
  geom_text(aes(x=10, y=2, label="amplitude", color="amplitude"),angle=90, size=6)+
  scale_color_gdocs()+
  theme_void()+
  theme(legend.position = "none")

ggsave("plots/ciclo_en.png")
```

tendencia poínomica

```{r eval=FALSE, include=TRUE}
nn = 200
#cambio_nivel
tendencia = exp(c(1:nn))
df <- data_frame(periodo=1:nn,tendencia=tendencia)

ggplot(df, aes(x=periodo, y=tendencia))+geom_line()

analyze.wavelet(my.data = df,
                my.series = "tendencia",
                verbose = F) %>%
  wt.image(.,color.key = "q",periodlab = "Largo del ciclo",timelab = "Tiempo calendario",
           plot.legend = FALSE,
           graphics.reset = FALSE)
```


Ahora grafico los wavelets de cada componente y la composición.

```{r message=FALSE, warning=FALSE, eval = FALSE}
analyze_and_save <- function(var, save = T, lang='es' ){
  
stopifnot(lang %in% c('es','en'))
  if (lang=='es') {
    periodlab = "Largo del ciclo"
    timelab = "Tiempo calendario"
    filename=glue("plots/espectograma_teorico_{var}.png")
  }
  if (lang=='en') {
    periodlab = "Length of the cycle"
    timelab = "Calendar time"
    filename=glue("plots/espectograma_teorico_{var}_en.png")
  }
  if (save) {
    png(filename, width = 1600, height = 1000, units = "p", pointsize = 40)
  }
  analyze.wavelet(my.data = df,
                my.series = var,
                verbose = F) %>%
  wt.image(.,color.key = "q",periodlab = periodlab ,timelab = timelab,
           plot.legend = FALSE,
           graphics.reset = FALSE)
  if (save) {
    dev.off()
  }
}

df <- data_frame(periodo=1:nn,impulso=impulso,tendencia=x0, ciclo_3=x3,ciclo_10=x10, ciclo_50=x50,ruido, composicion_series=x)

variables <- c("impulso","tendencia", "ciclo_3","ciclo_10", "ciclo_50","ruido", "composicion_series")
for (var in variables) {
  analyze_and_save(var)
}

for (var in variables) {
  analyze_and_save(var, lang = 'en')
}

```

- _cte_, que representa un cambio de nivel en 500, muestra todas las frecuencias de onda sobre la vertical en el punto de cambio de nivel. 
- La tendencia no se puede representar bien en el espectograma, porque el comportamiento ciclico es nulo. La diferencia con el ruido es que pasa de valores de intensidad más bajos al principio del período y tiempo (los ejes) a valores más altos. El ruido normal muestra mayores valores en las frecuencias mas altas (los periodos más bajos).
- Cada uno de los 3 componentes cíclicos se representa a la altura  de su period ()
- La amplitud del ciclo se representa con la intensidad (representado en la escala cromática)


La resolución del espectograma depende de la cantidad de observaciones

```{r message=FALSE, warning=FALSE}
nn = 25
cte= rep(100,nn)
x0 = ts(c(1:nn)/2)
x3 = ts(10*sin((2*pi/3)*c(1:nn)))
x10 = ts(20*sin((2*pi/10)*c(1:nn)))
x50 = ts(30*sin((2*pi/50)*c(1:nn)))
x = cte+ x0+  x3 + x10 + x50

df <- data_frame(period=1:nn,cte=cte,x0=x0, x3=x3,x10=x10, x50=x50, x=x)


```

- En x10 se ve que cuando el ciclo es demasiado grande respecto a la cantidad de observaciones, se satura en la parte superior de _period_


## Wavelet Base

```{r}
# install.packages("Rwave")
ejemplo_1 = Rwave::morlet(300,150,10,w0 = 2*pi) %>% Re(.)
ejemplo_2 = Rwave::morlet(300,150,20,w0 = 2*pi) %>% Re(.)
ejemplo_3 = Rwave::morlet(300,150,50,w0 = 2*pi) %>% Re(.)

ejemplo_escala =  data.frame(ejemplo_1,ejemplo_2,ejemplo_3, time =1:length(ejemplo_1)) %>%
  gather(wavelet,value,1:3) %>% 
  ggplot(.,aes(time,value, color=wavelet))+
  geom_line(size=1)+
  theme_minimal()+
    facet_grid(.~wavelet)+
  labs(title= "Transformaciones de escala", x="Tiempo", y="")+
  theme(legend.position = "none",
        strip.text = element_blank())

ejemplo_1 = Rwave::morlet(300,50,15,w0 = 2*pi) %>% Re(.)
ejemplo_2 = Rwave::morlet(300,150,15,w0 = 2*pi) %>% Re(.)
ejemplo_3 = Rwave::morlet(300,250,15,w0 = 2*pi) %>% Re(.)

ejemplo_traslacion= data.frame(ejemplo_1,ejemplo_2,ejemplo_3, time =1:length(ejemplo_1)) %>%
  gather(wavelet,value,1:3) %>% 
  ggplot(.,aes(time,value, color=wavelet))+
  geom_line(size=1)+
  theme_minimal()+
  labs(title= "Traslaciones", x="Tiempo", y="")+
      facet_grid(wavelet~.)+
 theme(legend.position = "none",
        strip.text = element_blank())

plot_grid(ejemplo_traslacion,ejemplo_escala,ncol=1)

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

```


```{r}

ejemplo_1 = Rwave::morlet(300,150,10,w0 = 2*pi) %>% Re(.)
ejemplo_2 = Rwave::morlet(300,150,20,w0 = 2*pi) %>% Re(.)
ejemplo_3 = Rwave::morlet(300,150,50,w0 = 2*pi) %>% Re(.)

ejemplo_escala =  data.frame(ejemplo_1,ejemplo_2,ejemplo_3, time =1:length(ejemplo_1)) %>%
  gather(wavelet,value,1:3) %>% 
  ggplot(.,aes(time,value, color=wavelet))+
  geom_line(size=1)+
  theme_minimal()+
    facet_grid(.~wavelet)+
  labs(title= "Scale transformations", x="Time", y="")+
  theme(legend.position = "none",
        strip.text = element_blank())


ejemplo_1 = Rwave::morlet(300,50,15,w0 = 2*pi) %>% Re(.)
ejemplo_2 = Rwave::morlet(300,150,15,w0 = 2*pi) %>% Re(.)
ejemplo_3 = Rwave::morlet(300,250,15,w0 = 2*pi) %>% Re(.)

ejemplo_traslacion= data.frame(ejemplo_1,ejemplo_2,ejemplo_3, time =1:length(ejemplo_1)) %>%
  gather(wavelet,value,1:3) %>% 
  ggplot(.,aes(time,value, color=wavelet))+
  geom_line(size=1)+
  theme_minimal()+
  labs(title= "Translations", x="Time", y="")+
      facet_grid(wavelet~.)+
 theme(legend.position = "none",
        strip.text = element_blank())

plot_grid(ejemplo_traslacion,ejemplo_escala,ncol=1)

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


# Series reales

### pbi y W en oro

```{r}

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)`,
         date = 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)`,
            date = parse_date_time(Year,"y")) %>%
  na.omit()

```

### salario

```{r message=FALSE, warning=FALSE}
wage_wavelet <- analyze.wavelet(my.data = wg_in_gold,
                                my.series = "value",
                                date.format = "%Y-%m-%d",
                                date.tz = "UTC", verbose = F)

```

```{r}
# png("plots/espectograma_wg.png", width = 1600, height = 1000, units = "p", pointsize = 40)
wt.image(wage_wavelet,color.key = "q", show.date = TRUE,
         label.time.axis = TRUE,
         plot.legend = FALSE,
         main = "Espectograma salario en oro")

# dev.off()
```

```{r}
png("plots/espectograma_wg_en.png", width = 1600, height = 1000, units = "p", pointsize = 40)
wt.image(wage_wavelet,color.key = "q", show.date = TRUE,
         label.time.axis = TRUE,
         plot.legend = FALSE)

dev.off()
```


log wage


```{r message=FALSE, warning=FALSE}
wage_wavelet <- wg_in_gold %>% 
  mutate(logvalue = log(value,base = 10)) %>% 
  analyze.wavelet(my.data = .,
                  my.series = "logvalue",
                  date.format = "%Y-%m-%d",
                  date.tz = "UTC", verbose = F)


# png("plots/espectograma_log_wg.png", width = 1600, height = 1000, units = "p", pointsize = 40)
wt.image(wage_wavelet,color.key = "q", show.date = TRUE,
         label.time.axis = TRUE,
         plot.legend = FALSE,
         main = "Espectograma salario en log(oro)")

# dev.off()

```


```{r}
png("plots/espectograma_log_wg_en.png", width = 1600, height = 1000, units = "p", pointsize = 40)

wt.image(wage_wavelet,color.key = "q", show.date = TRUE,
         label.time.axis = TRUE,
         plot.legend = FALSE)

dev.off()
```

### PBI

```{r}
gdp_vars <- gdp_in_gold %>% 
  mutate(var = ((value-lag(value))/lag(value))) %>% 
  filter(!is.na(var))


gdp_wavelet <- analyze.wavelet(my.data = gdp_vars,
                                my.series = "var",
                                date.format = "%Y-%m-%d",
                                date.tz = "UTC")
wt.image(gdp_wavelet,color.key = "q",
         show.date = TRUE,
         plot.legend = FALSE,
         label.time.axis = TRUE,
         main = "Espectograma vars PBI en oro")
```


```{r}

gdp_wavelet <- analyze.wavelet(my.data = gdp_in_gold,
                                my.series = "value",
                                date.format = "%Y-%m-%d",
                                date.tz = "UTC")

# png("plots/espectograma_gdp.png", width = 1600, height = 1000, units = "p", pointsize = 40)
wt.image(gdp_wavelet,color.key = "q",
         show.date = TRUE,
         plot.legend = FALSE,
         label.time.axis = TRUE,
         main = "Espectograma PBI en oro")
# dev.off()

```

```{r}
png("plots/espectograma_gdp_en.png", width = 1600, height = 1000, units = "p", pointsize = 40)
wt.image(gdp_wavelet,color.key = "q",
         show.date = TRUE,
         plot.legend = FALSE,
         label.time.axis = TRUE)
dev.off()

```

#### Escala log


```{r}

gdp_wavelet <- gdp_in_gold %>%
  mutate(logvalue = log(value,base = 10)) %>% 
  analyze.wavelet(my.data = .,
                  my.series = "logvalue",
                  date.format = "%Y-%m-%d",
                  date.tz = "UTC")  

# png("plots/espectograma_log_gdp.png", width = 1600, height = 1000, units = "p", pointsize = 40)

wt.image(gdp_wavelet,color.key = "q",
         periodlab="Frecuencia de ciclo",
         label.period.axis =TRUE,
         plot.legend = FALSE,
         spec.period.axis = list(at = c(3,7,20,50), labels = TRUE),
         n.levels = 100,
         show.date = TRUE,
         label.time.axis = TRUE, 
         main = "Espectograma log(PBI en oro)")
# dev.off()

```

```{r}
png("plots/espectograma_log_gdp_en.png", width = 1600, height = 1000, units = "p", pointsize = 40)

wt.image(gdp_wavelet,color.key = "q",
         label.period.axis =TRUE,
         plot.legend = FALSE,
         spec.period.axis = list(at = c(3,7,20,50), labels = TRUE),
         n.levels = 100,
         show.date = TRUE,
         label.time.axis = TRUE)
dev.off()

```


observaciones:

- En escala logarítima se reduce la heterocedasticidad de la serie, y sue puede observar los ciclos en un período más extendido del tiempo. especialmente el de 50 años
- Más allá de que hasta el 1900 la serie no tenga demasiada información, si filtarmos para quedarnos sólo con el último siglo, se perdería resolución y no se podría ver el ciclo de 50 años.


## otras series

### Tasa de interés de largo plazo

```{r}
ir <- ir %>% 
  mutate(date = parse_date_time(Year,"y")) %>% 
  select(date, value = `US Long-Term Rate: Consistent Series`)

ir_wavelet <- analyze.wavelet(my.data = ir,
                                my.series = "value",
                                date.format = "%Y-%m-%d",
                                date.tz = "UTC", verbose = F)

# png("plots/espectograma_ir.png", width = 1600, height = 1000, units = "p", pointsize = 40)
wt.image(ir_wavelet,color.key = "q", show.date = TRUE,
         label.time.axis = TRUE,
         plot.legend = FALSE,
         main = "Espectograma Tasa de interés de largo plazo")

# dev.off()

```

### IPC

```{r}
ipc <- ipc %>% 
  mutate(date = parse_date_time(Year,"y")) %>% 
  select(date, value = `U.S. Consumer Price Index *`)
ipc_wavelet <- analyze.wavelet(my.data = ipc,
                                my.series = "value",
                                date.format = "%Y-%m-%d",
                                date.tz = "UTC", verbose = F)

# png("plots/espectograma_ipc.png", width = 1600, height = 1000, units = "p", pointsize = 40)
wt.image(ipc_wavelet,color.key = "q", show.date = TRUE,
         label.time.axis = TRUE,
         plot.legend = FALSE,
         main = "Espectograma indice de precios al consumidor")

# dev.off()
```

### PBI real p/c

```{r}
real_gdp <- gdp %>% 
  mutate(date = parse_date_time(Year,"y")) %>% 
  select(date, value = `Real GDP per capita (year 2012 dollars)`)
real_gdp_wavelet <- analyze.wavelet(my.data = real_gdp,
                                my.series = "value",
                                date.format = "%Y-%m-%d",
                                date.tz = "UTC", verbose = F)

# png("plots/espectograma_real_gdp.png", width = 1600, height = 1000, units = "p", pointsize = 40)
wt.image(ipc_wavelet,color.key = "q", show.date = TRUE,
         label.time.axis = TRUE,
         plot.legend = FALSE,
         main = "Espectograma PBI real per cápita")

# dev.off()
```


#### pib real
```{r}
real_gdp <- gdp %>% 
  mutate(date = parse_date_time(Year,"y")) %>% 
  select(date, value = `Real GDP (millions of 2012 dollars)`)
real_gdp_wavelet <- analyze.wavelet(my.data = real_gdp,
                                my.series = "value",
                                date.format = "%Y-%m-%d",
                                date.tz = "UTC", verbose = F)



wt.image(ipc_wavelet,color.key = "q", show.date = TRUE,
         label.time.axis = TRUE,
         plot.legend = FALSE,
         main = "Espectograma PBI real")



```



### UK

```{r}
uk_gdp <- uk_gdp %>% 
  mutate(date = parse_date_time(Year,"y"),
         log_gdp = log(gdp_in_gold)) 


uk_gdp_filt <- uk_gdp %>% 
  filter(Year %in% c(1700:1900))

gdp_wavelet <- analyze.wavelet(my.data = uk_gdp_filt,
                                my.series = "gdp_in_gold",
                                date.format = "%Y-%m-%d",
                                date.tz = "UTC")

# png("plots/espectograma_gdp_uk.png", width = 1600, height = 1000, units = "p", pointsize = 40)
wt.image(gdp_wavelet,color.key = "q",
         show.date = TRUE,
         plot.legend = FALSE,
         label.time.axis = TRUE,
         main = "Espectograma PBI UK en oro")
# dev.off()
```

```{r}
png("plots/espectograma_gdp_uk_en.png", width = 1600, height = 1000, units = "p", pointsize = 40)
wt.image(gdp_wavelet,color.key = "q",
         show.date = TRUE,
         plot.legend = FALSE,
         label.time.axis = TRUE)
dev.off()
```




```{r}
gdp_wavelet <- analyze.wavelet(my.data = uk_gdp_filt,
                                my.series = "log_gdp",
                                date.format = "%Y-%m-%d",
                                date.tz = "UTC")

# png("plots/espectograma_log_gdp_uk.png", width = 1600, height = 1000, units = "p", pointsize = 40)
wt.image(gdp_wavelet,color.key = "q",
         show.date = TRUE,
         plot.legend = FALSE,
         label.time.axis = TRUE,
         main = "Espectograma log(PBI UK en oro)")
# dev.off()
```

```{r}
png("plots/espectograma_log_gdp_uk_en.png", width = 1600, height = 1000, units = "p", pointsize = 40)
wt.image(gdp_wavelet,color.key = "q",
         show.date = TRUE,
         plot.legend = FALSE,
         label.time.axis = TRUE)
dev.off()
```


serie centrada

```{r}

gdp_var <- uk_gdp %>% 
  select(Year, gdp_in_gold,date) %>%
  # filter(Year<1950) %>% 
  mutate(variacion = ((gdp_in_gold - lag(gdp_in_gold))/lag(gdp_in_gold))) %>%
  filter(!is.na(variacion))


gdp_wavelet <- analyze.wavelet(my.data = gdp_var,
                                my.series = "variacion",
                                date.format = "%Y-%m-%d",
                                date.tz = "UTC")

# png("plots/espectograma_gdp.png", width = 1600, height = 1000, units = "p", pointsize = 40)
wt.image(gdp_wavelet,color.key = "q",
         show.date = TRUE,
         plot.legend = FALSE,
         label.time.axis = TRUE,
         main = "Espectograma variaciones PBI en oro")
# dev.off()


```


### tendencia

```{r}

uk_gdp_filt$index <- 1:nrow(uk_gdp_filt)
loess_gdp <- loess(gdp_in_gold~index, data = uk_gdp_filt, span = 0.75)

uk_gdp_filt <- uk_gdp_filt %>% 
  mutate(trend = predict(loess_gdp),
         gdp_detrend = gdp_in_gold -trend)

ggplot(uk_gdp_filt, aes(x=date ))+
  geom_line(aes(y= gdp_in_gold, color = "serie original"), size=1)+
  geom_line(aes(y=trend, color = "tendencia loess"), size=1)+
  geom_line(aes(y= gdp_detrend, color = "sin tendencia"), size=1)+
  # scale_x_datetime(date_breaks = "25 years",labels = date_format("%Y") )+  
    theme_minimal()+
  labs(x="", y="PBI en oro", title= "Producto Bruto Interno Reino Unido",
       subtitle= "Millones de onzas de oro, 1700-1900. Tendencia")+
  theme(text = element_text(size = 20),
        legend.position = "bottom",
        legend.title = element_blank())

ggsave("plots/pbi_uk_tendencias.png", scale=1)
```

```{r}

# uk_gdp$index<- 1:nrow(uk_gdp)
# loess_gdp <- loess(gdp_in_gold~index, data = uk_gdp, span = 0.75)
# uk_gdp %>% 
#   mutate(trend = predict(loess_gdp),
#          gdp_detrend = gdp_in_gold -trend) %>% 

# ggplot(., aes(x=date ))+
ggplot(uk_gdp_filt, aes(x=date ))+
  geom_line(aes(y= gdp_in_gold, color = "original series"), size=1)+
  geom_line(aes(y=trend, color = "loess trend"), size=1)+
  geom_line(aes(y= gdp_detrend, color = "detrend"), size=1)+
  theme_minimal()+
  labs(x="", y="GDP in gold")+
  theme(text = element_text(size = 20),
        legend.position = "bottom",
        legend.title = element_blank())

ggsave("plots/pbi_uk_tendencias_en.png", scale=1)
```

```{r}
# sin tendencia

gdp_wavelet <- analyze.wavelet(my.data = uk_gdp_filt,
                                my.series = "gdp_detrend",
                                date.format = "%Y-%m-%d",
                                date.tz = "UTC")

png("plots/espectograma_gdp_uk_sinTend.png", width = 1600, height = 1000, units = "p", pointsize = 40)

wt.image(gdp_wavelet,color.key = "q",
         show.date = TRUE,
         plot.legend = FALSE,
         label.time.axis = TRUE,
         main = "Espectograma PBI UK en oro, sin tendencia")
dev.off()

```


```{r}
png("plots/espectograma_gdp_uk_sinTend_en.png", width = 1600, height = 1000, units = "p", pointsize = 40)

wt.image(gdp_wavelet,color.key = "q",
         show.date = TRUE,
         plot.legend = FALSE,
         label.time.axis = TRUE)
dev.off()
```


```{r}
# tendencia
gdp_wavelet <- analyze.wavelet(my.data = uk_gdp_filt,
                                my.series = "trend",
                                date.format = "%Y-%m-%d",
                                date.tz = "UTC")


png("plots/espectograma_gdp_uk_Tend.png", width = 1600, height = 1000, units = "p", pointsize = 40)

wt.image(gdp_wavelet,color.key = "q",
         show.date = TRUE,
         plot.legend = FALSE,
         label.time.axis = TRUE,
         main = "Espectograma PBI UK en oro, tendencia")
dev.off()
```

```{r}
png("plots/espectograma_gdp_uk_Tend_en.png", width = 1600, height = 1000, units = "p", pointsize = 40)

wt.image(gdp_wavelet,color.key = "q",
         show.date = TRUE,
         plot.legend = FALSE,
         label.time.axis = TRUE)
dev.off()
```




### oro


```{r}
gold <- gold %>% 
  mutate(date = parse_date_time(Year, "y"),
         value= `New York Market Price (U.S. dollars per fine ounce)`)

gold_wavelet <- analyze.wavelet(my.data = gold,
                                my.series = "value",
                                date.format = "%Y-%m-%d",
                                date.tz = "UTC")

# png("plots/espectograma_gdp.png", width = 1600, height = 1000, units = "p", pointsize = 40)
wt.image(gold_wavelet,color.key = "q",
         show.date = TRUE,
         plot.legend = FALSE,
         label.time.axis = TRUE,
         main = "Espectograma oro")

```




## cross-wavelet

```{r}

df <- left_join(gold, ipc) %>% 
  select(Year, oro =`New York Market Price (U.S. dollars per fine ounce)`, 
           ipc = `U.S. Consumer Price Index *`) %>% 
  mutate(oro = ts(oro),
         ipc= ts(ipc),
         date = parse_date_time(Year, "y")) %>% 
  left_join(gdp_in_gold) %>%
  rename(gdp = value) %>% 
              # select(Year, gdp=`Nominal GDP (million of Dollars)`)) %>% 
  filter(!is.na(date), !is.na(gdp))

##IPC vs oro
df$index <- 1:nrow(df)

loess_oro <- loess(oro~index, data = df, span = 0.75)
loess_ipc <- loess(ipc~index, data = df, span = 0.75)
loess_gdp <- loess(gdp~index, data = df, span = 0.75)

df <- df %>% 
  mutate(trend_oro = predict(loess_oro),
         trend_ipc = predict(loess_ipc),
         trend_gdp = predict(loess_gdp),
         oro_detrend = oro - trend_oro,
         ipc_detrend = ipc - trend_ipc,
         gdp_detrend = gdp - trend_gdp)

```

### espectograma detrended

```{r}
gold_wavelet <- analyze.wavelet(my.data = df,
                                my.series = "gdp_detrend",
                                date.format = "%Y-%m-%d",
                                date.tz = "UTC")

# png("plots/espectograma_gdp.png", width = 1600, height = 1000, units = "p", pointsize = 40)
wt.image(gold_wavelet,color.key = "q",
         show.date = TRUE,
         plot.legend = FALSE,
         label.time.axis = TRUE,
         main = "Espectograma gdp_detrend")

```



cross wavelets

my.pair = c("oro_detrend", "ipc_detrend"),

```{r, eval = F}
my.wc <- analyze.coherency(df,
                           my.pair = c("oro_detrend", "ipc_detrend"),
                           loess.span =0,
                           date.format = "%Y-%m-%d",
                           date.tz = "UTC"
                           )
wc.image(my.wc,
         color.key = "q",
         plot.legend = FALSE,
         label.time.axis = TRUE,
         show.date = TRUE,
         main = "cross-wavelet IPC-oro")

```

c("gdp_detrend","ipc_detrend"),
```{r, eval = F}
my.wc <- analyze.coherency(df,
                           my.pair = c("gdp_detrend","ipc_detrend"),
                           loess.span =0,
                           date.format = "%Y-%m-%d",
                           date.tz = "UTC"
                           )

wc.image(my.wc,
         color.key = "q",
         plot.legend = FALSE,
         label.time.axis = TRUE,
         show.date = TRUE,
         main = "cross-wavelet PBI-IPC")
```
c("gdp_detrend","oro_detrend"),
```{r, eval = F}
my.wc <- analyze.coherency(df,
                           my.pair = c("gdp_detrend","oro_detrend"),
                           loess.span =0,
                           date.format = "%Y-%m-%d",
                           date.tz = "UTC"
                           )

wc.image(my.wc,
         color.key = "q",
         plot.legend = FALSE,
         label.time.axis = TRUE,
         show.date = TRUE,
         main = "cross-wavelet PBI-oro")
```