En el siguiente post se vera como podemos descargar datos directamente del BCRP usando el Web Scraping y posteriormente se usara los difrenetes algoritmos de regresión los cuales son:

Además se vera como corregir la autocorrelación en la parte final.

Web Scraping

Paquetes a utilizar

ipak <- function(pkg){
  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
  if (length(new.pkg)) 
    install.packages(new.pkg, dependencies = TRUE)
  sapply(pkg, require, character.only = TRUE)
}
packages <- c("tidyverse","data.table", "XML","httr","xts")
ipak(packages)

Función de extracción

extraccion<- function(codigo=NULL,start=NULL,end=NULL){
  vector<- vector()
  datos<- NULL
  for(i in 1:length(codigo)){
                               ex<-  paste0("https://estadisticas.bcrp.gob.pe/estadisticas/series/api/", 
                                                       codigo[i], "/xml/",start,"/",end)
                               
  
                               ObjXML <- httr::GET(ex)
                               L_parseXML <- xmlParse(ObjXML)
                              
                               L_XML <- xmlToList(L_parseXML)
                               
                               dt_BCRP <- data.table::data.table()
                               for (a in 1:length(L_XML$periods)){
                                                                   dt_temp <- data.table(
                                                                                         serie= L_XML$periods[[a]]$v)
                                                                   dt_BCRP <- rbind(dt_BCRP,dt_temp)
                               }
                               vector[i]<- dt_BCRP[,serie:=as.numeric(serie)] 
  }
  for(i in  1: length(codigo)){
    datos<- cbind(datos,vector[[i]])
  }
  datos<- as.data.frame(datos)
  names(datos)<- codigo
  return(datos)
  }

Prueba de la función

datos<-extraccion(c("PN01488BM","PN01210PM" ),"2000-1","2007-2")
names(datos)<- c("Exportaciones","Tipo_de_cambio")
attach(datos)

Represenrtación de los resultados

datos  %>%round(digits = 2) %>% DT::datatable()

Analísis de Regresión

Regresión Simple

datos=datos
modelo_r<- lm(Exportaciones~Tipo_de_cambio,data = datos)
linea<- predict(modelo_r,datos)    
a1<-ggplot(datos,mapping = aes(y=Exportaciones,x=Tipo_de_cambio))+
      geom_point(color="red")+geom_line(mapping = aes(x=Tipo_de_cambio,
                                             y=linea))+theme_minimal()+
   xlab("Tipo de cambio")+labs(title ="Regresión lineal")
E_c_M_1= (Exportaciones-linea)^2 %>% mean()
a1

Regresión polinomial

modelo_r<- lm(Exportaciones~poly(Tipo_de_cambio,2),data = datos)
    linea_1<- predict(modelo_r,datos)    
    a2<-ggplot(datos,mapping = aes(y=Exportaciones,x=Tipo_de_cambio))+
      geom_point(color="red")+geom_line(mapping = aes(x=Tipo_de_cambio,
                                                      y=linea_1))+theme_minimal()+
      xlab("Tipo de cambio")+labs(title ="Regresión polinomica")
    E_c_M_2= (Exportaciones-linea_1)^2 %>% mean()
a2

Maquinas de soporte vectorial

library(e1071)
    modelo_r<- svm(Exportaciones~Tipo_de_cambio,data = datos)
    linea_2<- predict(modelo_r,datos)    
    a3<-ggplot(datos,mapping = aes(y=Exportaciones,x=Tipo_de_cambio))+
      geom_point(color="red")+geom_line(mapping = aes(x=Tipo_de_cambio,
                                                      y=linea_2))+theme_minimal()+
    xlab("Tipo de cambio")+labs(title ="Maquina de soporte vectorial")
    
    E_c_M_3= (Exportaciones-linea_2)^2 %>% mean()
a3

Arboles de decisión

library(rpart)
    modelo_r<- rpart(Exportaciones~Tipo_de_cambio,data = datos)
    linea_3<- predict(modelo_r,datos)    
    a4<-ggplot(datos,mapping = aes(y=Exportaciones,x=Tipo_de_cambio))+
      geom_point(color="red")+geom_line(mapping = aes(x=Tipo_de_cambio,
                                                      y=linea_3))+theme_minimal()+
      xlab("Tipo de cambio")+labs(title ="Arboles de desición")
    E_c_M_4= mean((Exportaciones-linea_3)^2)
  a4

Bosques aleatorios

# library(randomForest)
    modelo_r<- randomForest(Exportaciones~Tipo_de_cambio,data = datos)
    linea_4<- predict(modelo_r,datos)    
    a5<-ggplot(datos,mapping = aes(y=Exportaciones,x=Tipo_de_cambio))+
      geom_point(color="red")+geom_line(mapping = aes(x=Tipo_de_cambio,
                                                      y=linea_4))+theme_minimal()+
      xlab("Tipo de cambio")+labs(title ="Bosques Aleatorios")
    E_c_M_5= (Exportaciones-linea_4)^2 %>% mean()    
a5

Elección del mejor modelo

Nota: El mejor modelo es en el cual se aplico el metódo de Random forest pues presente un menor \(R^2\).

Corrección de la autocorrelación

Para ver si presenta autocorrelación usamos el test de durwin watson que nos permite ver si se presenta esto. Para corregirla se presenta las siguientes funciones.

Uso de la prueba de autocorrelación

modelo_r<- lm(Exportaciones~Tipo_de_cambio,data = datos)
summary(modelo_r)
## 
## Call:
## lm(formula = Exportaciones ~ Tipo_de_cambio, data = datos)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -536.5 -200.2  -24.1  190.2 1006.0 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     14778.8      949.9   15.56   <2e-16 ***
## Tipo_de_cambio  -4021.5      277.7  -14.48   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 286.6 on 84 degrees of freedom
## Multiple R-squared:  0.714,  Adjusted R-squared:  0.7106 
## F-statistic: 209.7 on 1 and 84 DF,  p-value: < 2.2e-16
b<- modelo_r$residuals
acf(b)

pacf(b)

Resultado: podemos ver que de forma gráfica se presenta autocorrelación

Ahora aplicamos una prueba formal

lmtest::dwtest(modelo_r) 
## 
##  Durbin-Watson test
## 
## data:  modelo_r
## DW = 0.5182, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

Podemos ver que el DW calculado es mayor a lo permitido así que podemos deducir que si presenta autocorrelación

Corrección

Agregando un resago

library(dynlm)      
modelo_a<- dynlm(Exportaciones~Tipo_de_cambio+L(Exportaciones))      
summary(modelo_a)
## Warning in summary.lm(modelo_a): essentially perfect fit: summary may be
## unreliable
## 
## Time series regression with "numeric" data:
## Start = 1, End = 86
## 
## Call:
## dynlm(formula = Exportaciones ~ Tipo_de_cambio + L(Exportaciones))
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.073e-13 -2.086e-14 -9.510e-15 -6.000e-17  7.833e-13 
## 
## Coefficients:
##                    Estimate Std. Error    t value Pr(>|t|)    
## (Intercept)       6.417e-12  5.901e-13  1.087e+01   <2e-16 ***
## Tipo_de_cambio   -1.762e-12  1.637e-13 -1.076e+01   <2e-16 ***
## L(Exportaciones)  1.000e+00  3.440e-17  2.907e+16   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.038e-14 on 83 degrees of freedom
## Multiple R-squared:      1,  Adjusted R-squared:      1 
## F-statistic: 1.477e+33 on 2 and 83 DF,  p-value: < 2.2e-16
c<- modelo_a$residuals

Analizando los resultados:

acf(c)

pacf(c)

lmtest::dwtest(modelo_a) 
## 
##  Durbin-Watson test
## 
## data:  modelo_a
## DW = 2.0655, p-value = 0.5445
## alternative hypothesis: true autocorrelation is greater than 0

Nota: Ahora este indicador se aproxima a 2 así que podemos decir que no presenta acutocrrelación.

Algortimo de ourcurt

library(orcutt)      
## Loading required package: lmtest
cor<-cochrane.orcutt(modelo_r)
summary.orcutt(cor)
## Call:
## lm(formula = Exportaciones ~ Tipo_de_cambio, data = datos)
## 
##                Estimate Std. Error t value Pr(>|t|)
## (Intercept)     1021.77    2150.97   0.475   0.6360
## Tipo_de_cambio   137.68     629.46   0.219   0.8274
## 
## Residual standard error: 165.273 on 83 degrees of freedom
## Multiple R-squared:  6e-04 ,  Adjusted R-squared:  -0.0115
## F-statistic: 0 on 1 and 83 DF,  p-value: < 8.274e-01
## 
## Durbin-Watson statistic 
## (original):    0.51820 , p-value: 1.866e-17
## (transformed): 2.77543 , p-value: 9.999e-01