TALLER 6 - ECONOMETRÍA 2

DESARROLLO DEL TALLER 6 DE ECONOMETRÍA 2

A continuación se muestran los códigos necesarios para la estimación de los modelos del taller.

Código del Desarrollo del Punto 3

#|----------------------------------------------------------------------------|#
#|-------->>| METODOLOGIA CLÁSICA -- SUAVIZAMIENTO - Holt-Winters |<<---------|#
#|----------------------------------------------------------------------------|#

#|-------------------------->>| CARGA DE PAQUETES |<<-------------------------|#

library(xts)
library(tseries)
library(forecast)
library(lmtest)
library(readxl)
library(graphics)
library(ggplot2)
library(stargazer)
library(writexl)
library(WriteXLS)


#|----------------------------------------------------------------------------|#

#| Directorio
setwd("C:/Programacion en R/5. ECONOMETRIA 2/Taller-6")

#| Datos como CT
DATA <- read_excel("DATA3.xlsx")
View(DATA)

#| Conversión a ST
y <- ts(DATA[,2], start=c(1956,1), frequency = 4)
y

#|----------------------------------------------------------------------------|#

#| Gráfica de la ST
par(mfrow=c(1,1))
plot(y, main="Producción de Cerveza en Australia", 
     col  = "black", 
     lwd  = 2,
     xlab = "Trimestres",
     ylab = "Producción de Cerveza")

#| Caracteristicas de la ST
class(y)          #| ts -> time serie
frequency(y)      #| 4


#|----------------------------------------------------------------------------|#

#| Efecto Estacional
par(mfrow=c(1,1))
boxplot(y ~ cycle(y),main="Producción de Cerveza en Australia", 
        lwd  = 2,
        xlab = "Trimestres",
        ylab = "Producción de Cerveza")

ggseasonplot(y, col=rainbow(12), year.labels=TRUE)+
  ylab("Producción de Cerveza") +
  xlab("Trimestre") +
  ggtitle("Producción de Cerveza en Australia cada Año por Trimestre")

ggseasonplot(y, year.labels=TRUE, continuous=TRUE)+
  ylab("Producción de Cerveza") +
  ggtitle("Gráfico de Evaluación de Estacionalidad  (1956-Q1 - 1994-Q2)")

#|----------------------->>| METODOLOGÍA BOX-JENKINS |<<----------------------|#

#| ESTIMACIÓN DE MODELOS CANDIDATOS A PGD SARIMA

#| SARIMA(0,1,2)(0,1,1)_4
pdqParamR = c(0, 1, 2)
pdqParamS = c(0, 1, 1)
MODELO1 <- arima(y, pdqParamR, seasonal = list(order = pdqParamS, period = 4))
MODELO1
coeftest(MODELO1) 

#| SARIMA(2,1,0)(1,1,0)_4
pdqParamR = c(2, 1, 0)
pdqParamS = c(1, 1, 0)
MODELO2 <- arima(y, pdqParamR, seasonal = list(order = pdqParamS, period = 4))
MODELO2
coeftest(MODELO2) 

#|----------------------------------------------------------------------------|#

#| DIAGNÓSTICO

#| Supuesto de NO AC
tsdiag(MODELO1) #| Modelo con NO AC
tsdiag(MODELO2) #| Modelo con AC

#| Supuesto de NORMALIDAD
jarque.bera.test(MODELO1$residuals) #| NORMAL
jarque.bera.test(MODELO2$residuals) #| NO NORMAL

#| Prueba de raíces inversas
autoplot(MODELO1) #| Estacionario e Invertible
autoplot(MODELO2) #| Estacionario e Invertible

#| Criterio Akaike y Bayesiano 

AIC(MODELO1)  #| Mejor Modelo 1 
AIC(MODELO2)

BIC(MODELO1)  #| Mejor Modelo 1
BIC(MODELO2)

#|----------------------------------------------------------------------------|#

#| Se pronostican 2 años, 24 meses por ser la serie mensual
par(mfrow=c(1,1))
MODELO1_PRED <- forecast(MODELO1, h=10)
plot(MODELO1_PRED,main="Pronóstico de la Producción de Cerveza en Australia", 
     lwd  = 2,
     xlab = "Trimestres",
     ylab = "Producción de Cerveza")

MODELO1_PRED

stargazer(MODELO1, type="text", out="ResultsPunto3.doc") 

ResultsPunto3a <- as.data.frame(MODELO1_PRED)
View(ResultsPunto3a)
write_xlsx(ResultsPunto3a,"ResultsPunto3a.xlsx")

#|----------------------->>| METODOLOGÍA HOLT-WINTERS|<<----------------------|#

#| HW estacional aditivo
hwa_y <- HoltWinters(y)
hwa_y
plot(hwa_y, col = 'black')

#| HW estacional multiplicativo
hwm_y <- HoltWinters(y, seasonal = "multiplicative")
hwm_y
plot(hwm_y, col = 'blue')



SCE_hwa_y <- hwa_y$SSE #| Suma Cuadrada de los Errores
ECM_hwa_y <- sqrt(hwa_y$SSE/154) #| Raíz del Error Cuadratico Medio

SCE_hwm_y <- hwm_y$SSE #| Suma Cuadrada de los Errores
ECM_hwm_y <- sqrt(hwm_y$SSE/154) #| Raíz del Error Cuadratico Medio

SCE_hwa_y 
SCE_hwm_y #| Mejor Modelo

ECM_hwa_y 
ECM_hwm_y 

#|----------------------------------------------------------------------------|#

#| Pronóstico
hwm_yf <- forecast(hwm_y, h=10, level = c(90,95), fan = TRUE)
hwm_yf
plot(hwm_yf,main="Pronóstico de la Producción de Cerveza en Australia", 
     lwd  = 2,
     xlab = "Trimestres",
     ylab = "Producción de Cerveza")

ResultsPunto3b <- as.data.frame(hwm_yf)
View(ResultsPunto3b)
write_xlsx(ResultsPunto3b,"ResultsPunto3b.xlsx")


#|----------------------------------------------------------------------------|#

#| Comparación de ambas metologías

par(mfrow=c(2,1))
plot(MODELO1_PRED,main="Pronóstico de Box-Jenkins", 
     lwd  = 2,
     xlab = "Trimestres",
     ylab = "Producción de Cerveza",
     xlim=c(1990,1995))
plot(hwm_yf,main="Pronóstico de Holt-Winters", 
     lwd  = 2,
     xlab = "Trimestres",
     ylab = "Producción de Cerveza",
     xlim=c(1990,1995))

Código del Desarrollo del Punto 4

#|----------------------------------------------------------------------------|#
#|----------->>| METODOLOGIA CLÁSICA -- DESAGREGACIÓN TEMPORAL |<<------------|#
#|----------------------------------------------------------------------------|#

#|-------------------------->>| CARGA DE PAQUETES |<<-------------------------|#
library(xts)
library(tseries)
library(forecast)
library(lmtest)
library(readxl)
library(ggplot2)
library(stargazer)
library(writexl)
library(WriteXLS)
library(tempdisagg)
library(xlsx)
library(tsbox)

#|----------------------------------------------------------------------------|#

#| Directorio
setwd("C:/Programacion en R/5. ECONOMETRIA 2/Taller-6")

#| Datos como CT
DATAa <- read_excel("DATA4a.xlsx")
View(DATAa)

#| Conversión a ST
Y_Trimestral <- ts(DATAa[,2], start=c(2005,1), frequency = 4)
Y_Trimestral


plot(Y_Trimestral)

#| El objetivo es volver la serie mensual

#| Datos como CT
DATAb <- read_excel("DATA4b.xlsx", sheet = 1)
View(DATAb)

#| Conversión a ST
VI <- ts(DATAb[,2], start=c(2005), frequency = 12)
VI

plot(VI)



#|----------------------------------------------------------------------------|#

#| CONVERSIÓN DE LA SERIE TRIMESTRAL A MENSUAL

#| Sin Variable Indicadora
M.MODEL1 <- td(Y_Trimestral ~ 1,      
               to = "monthly", 
               method = "fast")
summary(M.MODEL1)

Y_Mensual_NO_VI <- predict(M.MODEL1)
Y_Mensual_NO_VI

#| Sobreposición de las gráficas - SIN ESCALA
GRAPH_1 <- ts_plot(ts_scale(ts_c(Y_Trimestral,Y_Mensual_NO_VI)),
                   title = "PIB Real Trimestral y PIB Real Desagregado Mensual",
                   subtitle = "Sin Variable Indicadora",
                   ylab = "Escala Ajustada - Sobreposición de Series ")

#|----------------------------------------------------------------------------|#

#| CONVERSIÓN DE LA SERIE TRIMESTRAL A MENSUAL

#| Con Variable Indicadora
M.MODEL2 <- td(Y_Trimestral ~ VI)
summary(M.MODEL2)

Y_Mensual_VI <- predict(M.MODEL2)
Y_Mensual_VI

GRAPH_2 <- ts_plot(ts_scale(ts_c(Y_Trimestral,Y_Mensual_VI)),
                   title = "PIB Real Trimestral y PIB Real Desagregado Mensual",
                   subtitle = "Con Variable Indicadora",
                   ylab = "Escala Ajustada - Sobreposición de Series ")



ResultsPunto4a <- as.data.frame(Y_Mensual_NO_VI)
View(ResultsPunto4a)

ResultsPunto4b <- as.data.frame(Y_Mensual_VI)
View(ResultsPunto4b)

write_xlsx(ResultsPunto4a,"ResultsPunto4b.xlsx")
write_xlsx(ResultsPunto4b,"ResultsPunto4b.xlsx")


GRAPH_3 <- ts_plot(ts_scale(ts_c(Y_Trimestral,Y_Mensual_NO_VI,Y_Mensual_VI)),
                   title = "PIB Real Trimestral y PIB Real Desagregado Mensual",
                   subtitle = "Sin Variable Indicadora y Con Variable Indicadora",
                   ylab = "Escala Ajustada - Sobreposición de Series ")

Código del Desarrollo del Punto 5

#|----------------------------------------------------------------------------|#
#|----------->>| METODOLOGIA CLÁSICA -- DESAGREGACIÓN TEMPORAL |<<------------|#
#|----------------------------------------------------------------------------|#

#|-------------------------->>| CARGA DE PAQUETES |<<-------------------------|#

library(writexl)
library(WriteXLS)
library(xts)
library(tseries)
library(forecast)
library(lmtest)
library(mFilter)
library(readxl)

#|----------------------------------------------------------------------------|#

#| Directorio
setwd("C:/Programacion en R/5. ECONOMETRIA 2/Taller-6")
dir()

#| Datos como CT
DATA <- read_excel("DATA5.xlsx")
View(DATA)

#| Conversión a ST
y <- ts(DATA[,2], start = c(1960), frequency = 1)

#| Características de la ST
class(y)
start(y)
end(y)

#|----------------------------------------------------------------------------|#

#| Gráfica de la ST
ts_plot(y,title = "Comportamiento del PIB Anual",
        subtitle = "Desde 1960 a 2020")

#|----------------------------------------------------------------------------|#

#| Logaritmo
ly <- log(y)
#| Logaritmo porque se puede interpretar el ciclo en porcentaje (%)
#| ya que el C=Y-T y en este caso la resta del logaritmo del (y) observado
#| con el (y) tendencial se puede ver como la resta de dos logaritmos

ts_plot(ly,title = "Comportamiento del Logaritmo del PIB Anual",
        subtitle = "Desde 1960 a 2020")

#|----------------------------------------------------------------------------|#

#| Filtro HP 

#| Frecuencia Anual - Lambda=100
filtro_100 <- hpfilter(ly, freq = 100, type ="lambda")
plot(filtro_100) #| Gráfico del Filtro HP


ly_trend_100 <- filtro_100$trend
ly_trend_100
ly_cycle_100 <- filtro_100$cycle
ly_cycle_100


#| Frecuencia Anual - Lambda=5
filtro_5 <- hpfilter(ly, freq = 5, type ="lambda")
plot(filtro_5) #| Gráfico del Filtro HP


ly_trend_5 <- filtro_5$trend
ly_trend_5
ly_cycle_5 <- filtro_5$cycle
ly_cycle_5

#| Gráfico del Ciclo y la Tendencia de la ST
par(mfrow=c(1,2))

plot(ly_trend_100, main="Log(Y) - Tendencia-Lambda (100)", col="red",lwd = 2)
lines(ly, col="blue",lwd = 2)

plot(ly_trend_5, main="Log(Y) - Tendencia-Lambda (5)", col="red",lwd = 2)
lines(ly, col="blue",lwd = 2)

plot(ly_cycle_100, main="Log(Y) - Ciclo-Lambda (100)", col="darkgreen",lwd = 2)
plot(ly_cycle_5, main="Log(Y) - Ciclo-Lambda (5)", col="darkgreen",lwd = 2)


dataframe = data.frame(ly,ly_trend_100, ly_cycle_100,ly_trend_5, ly_cycle_5)
View(dataframe)

write_xlsx(dataframe,"ResultsPunto5.xlsx")

Código del Desarrollo del Punto 6

#|----------------------------------------------------------------------------|#
#|----------->>| METODOLOGIA CLÁSICA -- DESAGREGACIÓN TEMPORAL |<<------------|#
#|----------------------------------------------------------------------------|#

#|-------------------------->>| CARGA DE PAQUETES |<<-------------------------|#

library(writexl)
library(WriteXLS)
library(xts)
library(tseries)
library(forecast)
library(lmtest)
library(mFilter)
library(readxl)
library(haven)
library(car)
library(stargazer)

#|----------------------------------------------------------------------------|#

#| Directorio
setwd("C:/Programacion en R/5. ECONOMETRIA 2/Taller-6")

#| Datos como CT
DATA <- read_excel("DATA6.xlsx")
View(DATA)

#| Conversión a ST
GR <- ts(DATA[,3], start = c(1980), frequency = 1)
GR
UN <- ts(DATA[,2], start = c(1980), frequency = 1)
UN

#|----------------------------------------------------------------------------|#

#| Gráfica de la ST
par(mfrow=c(1,1))
ts.plot(GR,main="Crecimiento del PIB real", col = "red")
ts.plot(UN,main="Tasa de Desempleo ", col = "red")

#|----------------------------------------------------------------------------|#

#| Filtro HP 

#| Frecuencia Anual - Lambda=100
HP_FILTER_GR <- hpfilter(GR, freq = 100, type ="lambda")
plot(HP_FILTER_GR) #| Gráfico del Filtro HP

GR_TREND <- HP_FILTER_GR$trend
GR_TREND

GR_CYCLE <- HP_FILTER_GR$cycle
GR_CYCLE

#| Frecuencia Anual - Lambda=100
HP_FILTER_UN <- hpfilter(UN, freq = 100, type ="lambda")
plot(HP_FILTER_UN) #| Gráfico del Filtro HP

UN_TREND <- HP_FILTER_UN$trend
UN_TREND

UN_CYCLE <- HP_FILTER_UN$cycle
UN_CYCLE


#| Gráfico del Ciclo y la Tendencia de la ST
par(mfrow=c(1,2))
plot(GR_TREND, main="Tendencia del Crecimiento del PIB real", col="red")
plot(GR_CYCLE, main="Ciclo del Crecimiento del PIB real", col="blue")

#| Gráfico del Ciclo y la Tendencia de la ST
par(mfrow=c(1,2))
plot(UN_TREND, main="Tendencia de la Tasa de Desempleo", col="red")
plot(UN_CYCLE, main="Ciclo de la Tasa de Desempleo", col="blue")

#|----------------------------------------------------------------------------|#

#| Comparación Interesante entre el ciclo del PIB y su tasa de crecimiento

par(mfrow=c(1,1))
plot(GR, main="Crevimiento del PIB real", col="red")
lines(GR_TREND, col="blue")

#|----------------------------------------------------------------------------|#

#| Exportar datos a Excel
dataframe = data.frame(GR,GR_TREND,GR_CYCLE, UN,UN_TREND,UN_CYCLE)
View(dataframe)

write_xlsx(dataframe,"ResultsPunto6.xlsx")


#|----------------------------------------------------------------------------|#

#| MODELO ECONOMÉTRICO

setwd("C:/Programacion en R/5. ECONOMETRIA 2/Taller-6")

#| Datos como CT
DATAb <- read_excel("DATA6b.xlsx")
View(DATAb)

CYC <- ts(DATAb[,3], start = c(1981), frequency = 1)
CYC
UND <- ts(DATAb[,2], start = c(1981), frequency = 1)
UND

OKUN_MODEL = lm(UND~CYC-1)
summary(OKUN_MODEL)

stargazer(OKUN_MODEL,type="text", out="ResultsPunto6.doc")

#| Errores del Modelo
error=residuals(OKUN_MODEL)

#| Pruebas de Normalidad
jarque.bera.test(error) #| NORMAL
shapiro.test(error)     #| NORMAL                                   

#| Pruebas de No Autocorrelación
ac <- acf(error, col="red", lwd=2, main="Correlograma - Errores")    

#| Pruebas de Homoscedasticidad
bpagan <- bptest(OKUN_MODEL, ~I(CYC))
bpagan               #| HM
ncvTest(OKUN_MODEL)  #| HM

stargazer(bpagan, out="ResultsPunto6bpagan.doc")