TRABAJO FINAL DE SERIES DE TIEMPO

library(TSA) library(ggplot2) library(ggfortify) library(forecast) library(fpp2) library(forecast) library(fma) library(expsmooth) library(seasonal) library(urca) library(seasonal) library(zoo) library(dynlm) library(urca) library(foreing) library(lmtest) library(vars)

TABLA DE CONTENIDOS

1. Breve referencia a los datos, frecuencia, unidad de medida, fuente, etc

2. Graficar la serie de tiempo, analizar patrones y observaciones atípicas. Apoye su análisis con una descomposición clásica.

3. Si es necesario, utilizar transformación Box-Cox / logaritmos para estabilizar la varianza. Presentar grafica

4. Modelar la serie de tiempo elegida como función de una tendencia, tendencia cuadrática y/o medias estacionales. Presente el diagnóstico del modelo elegido.

5. En caso de estacionalidad de la serie de tiempo, aplicar diferencias estacionales. Presentar grafica

6. Usar prueba Dickey-Fuller para evaluar el orden de integración de la serie. Diferenciar hasta que la serie sea estacionaria.

7. Usar herramientas ACF,PACF,EACF y/o criterios de Akaike/Bayes para construir propuestas de modelos.

8. Como parte del diagnóstico del modelo, analizar los residuos y aplicar la prueba Ljung-Box.

9. Usar la función auto.arima() y compare sus resultados con el inciso anterior.

10. Presente la ecuación final y describa.

11. Realizar el pronóstico ARIMA para dos años.

12. Utilizando métodos de pronósticos simples y/o suavizamiento exponencial generar un pronóstico para dos años. Elegir el método que presente la raíz del error medio al cuadrado más pequeño.

13. Usando el an´alisis de transferencia, modelar observaciones cambios estructurales y/o observaciones atípicas en el contexto ARIMAX (punto extra, opcional)

14. Analizar un VAR, necesita incluir los puntos siguientes

a) Buscar otra variable (incluso m´as variables) y redactar la teoría que se encuentra detrás de la relación entre variables que propone.

b) Evaluar la estacionalidad de la nueva variable.

c) Proponga el orden de rezagos óptimo para el VAR y presentar la estimaci´on.

d) Probar la correlación serial y normalidad en los residuos.

e) Resultados de la prueba de causalidad de Granger, y

f) Crear un pronóstico para dos años usando el VAR especificado.

g) Interpretar las funciones impulso-respuesta.

15. Evaluar la potencial cointegraci´on usando la metodologíia de Engle-Granger y Johansen.

1. Breve referencia a los datos, frecuencia, unidad de medida, fuente, etc.

La tasa de desempleo, también conocida como tasa de paro, mide el nivel de desocupación en relación a la población activa. En otras palabras, es la parte de la población que estando en edad, condiciones y disposición de trabajar -población activa- no tiene puesto de trabajo.

La tasa de desempleo es muy útil para conocer las personas que no están trabajando. Su fórmula de cálculo es la población de 16 años y más que no está trabajando y busca trabajo, dividido entre la población económicamente activa de 16 años y más, esto es, ocupados más desocupados. Se calcula de la siguiente manera: tasa de desempleo=(No. de desempleados/Poblacion activa)*100 datos sacdos de http://www.fao.org

2. Graficar la serie de tiempo, analizar patrones y observaciones atípicas. Apoye su análisis con una descomposición clásica.

base<-read.csv(file.choose())
tdd<-ts(base$tdd, start= c(1987,1),frequency= 12)
autoplot(tdd)

ggseasonplot(tdd)

ggseasonplot(tdd, polar=TRUE)

aquí podemos apreciar como tenemos una deserie donde aument eldesempleo en los 90’s por el efecto tequila, donde tuvimos una depresicación de la moneda mexicana, de igual manera podemos observar como baja poco a poco hasta el año 2000 donde tenemos un incremento hasta llegar a otra crisis en el años 2008-2009 pero no tan fuerte como en los 90’s

fit <- decompose(tdd, type='multiplicative')
autoplot(fit)

Para el caso clásico obtenemos los componentes separados de nuestra serie, en relacion a la estacionalidad, tendencia y residuos.

fit2 <- seas(tdd, x11="")
autoplot(fit2) + ggtitle("X11 decomposicion de la tasa de desempleo")

Se realizo una gráfica que compara el comportamiento de nuestros datos con el de un ajuste estacional y observamos un comportamiento más preciso en el componente tendencial, además de que la estacionalidad refleja cambios que van acortando los extremos verticales con una pequeña alza a mitad del periodo, que después regresa a su actividad inicial.

3. Si es necesario, utilizar transformación Box-Cox / logaritmos para estabilizar la varianza. Presentar grafica

En nuestros caso se estan trabajando con índices, por lo cual no es correcto el hacer uso de logaritmos, por lo que solo consideramos la aplicación de la diferenciación. Por lo anterior dicho, no hacemos consideración de la transformación Box Cox.

4. Modelar la serie de tiempo elegida como función de una tendencia, tendencia cuadrática y/o medias estacionales. Presente el diagnóstico del modelo elegido.

modelo1 <- lm(tdd~time(tdd))
summary(modelo1)

Call:
lm(formula = tdd ~ time(tdd))

Residuals:
    Min      1Q  Median      3Q     Max 
-1.6330 -0.6457 -0.1770  0.5039  4.0496 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -77.438827  10.736534  -7.213 3.07e-12 ***
time(tdd)     0.040584   0.005361   7.570 2.94e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9403 on 374 degrees of freedom
Multiple R-squared:  0.1329,    Adjusted R-squared:  0.1305 
F-statistic: 57.31 on 1 and 374 DF,  p-value: 2.941e-13

ggplot(base, aes(y=tdd, x=time(tdd))) + geom_point() + geom_smooth(method = 'lm')

autoplot(modelo1, which = 1:6, ncol = 2, label.size = 3)

shapiro.test(rstudent(modelo1))

    Shapiro-Wilk normality test

data:  rstudent(modelo1)
W = 0.92678, p-value = 1.326e-12

autoplot(acf(rstudent(modelo1),plot=F))

mes. <- season(tdd)
modelo2 <- lm(tdd~mes.)
summary(modelo2)

Call:
lm(formula = tdd ~ mes.)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.7656 -0.7581 -0.2223  0.7524  3.7032 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)    3.818750   0.180821  21.119   <2e-16 ***
mes.February   0.046875   0.255720   0.183    0.855    
mes.March      0.015625   0.255720   0.061    0.951    
mes.April      0.031250   0.255720   0.122    0.903    
mes.May        0.039315   0.257774   0.153    0.879    
mes.June       0.039315   0.257774   0.153    0.879    
mes.July       0.048992   0.257774   0.190    0.849    
mes.August     0.078024   0.257774   0.303    0.762    
mes.September -0.015524   0.257774  -0.060    0.952    
mes.October   -0.018750   0.257774  -0.073    0.942    
mes.November   0.007056   0.257774   0.027    0.978    
mes.December  -0.063911   0.257774  -0.248    0.804    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.023 on 364 degrees of freedom
Multiple R-squared:  0.001313,  Adjusted R-squared:  -0.02887 
F-statistic: 0.04351 on 11 and 364 DF,  p-value: 1

ggplot(base, aes(y=tdd, x=time(tdd))) + geom_point() + geom_smooth(method = 'lm')

autoplot(modelo2, which = 1:6, ncol = 2, label.size = 3)

shapiro.test(rstudent(modelo2))

    Shapiro-Wilk normality test

data:  rstudent(modelo2)
W = 0.94926, p-value = 4.545e-10

autoplot(acf(rstudent(modelo2),plot=F))

modelo3 <- lm(tdd~mes. -1)
summary(modelo3)

Call:
lm(formula = tdd ~ mes. - 1)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.7656 -0.7581 -0.2223  0.7524  3.7032 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
mes.January     3.8187     0.1808   21.12   <2e-16 ***
mes.February    3.8656     0.1808   21.38   <2e-16 ***
mes.March       3.8344     0.1808   21.20   <2e-16 ***
mes.April       3.8500     0.1808   21.29   <2e-16 ***
mes.May         3.8581     0.1837   21.00   <2e-16 ***
mes.June        3.8581     0.1837   21.00   <2e-16 ***
mes.July        3.8677     0.1837   21.05   <2e-16 ***
mes.August      3.8968     0.1837   21.21   <2e-16 ***
mes.September   3.8032     0.1837   20.70   <2e-16 ***
mes.October     3.8000     0.1837   20.68   <2e-16 ***
mes.November    3.8258     0.1837   20.82   <2e-16 ***
mes.December    3.7548     0.1837   20.44   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.023 on 364 degrees of freedom
Multiple R-squared:  0.9356,    Adjusted R-squared:  0.9335 
F-statistic: 440.7 on 12 and 364 DF,  p-value: < 2.2e-16

ggplot(base, aes(y=tdd, x=time(tdd))) + geom_point() + geom_smooth(method = 'lm')

autoplot(modelo3, which = 1:6, ncol = 2, label.size = 3)

shapiro.test(rstudent(modelo3))

    Shapiro-Wilk normality test

data:  rstudent(modelo3)
W = 0.94926, p-value = 4.545e-10

autoplot(acf(rstudent(modelo3),plot=F))

5. En caso de estacionalidad de la serie de tiempo, aplicar diferencias estacionales. Presentar grafica

tdd.d<-autoplot(diff(tdd))

autoplot(tdd.d)+ ggtitle("tasa de desempleo")

6. Usar prueba Dickey-Fuller para evaluar el orden de integración de la serie. Diferenciar hasta que la serie sea estacionaria.

Diferenciar hasta que la serie sea estacionaria. Ho: Existe raíz unitaria Ha: No existe raíz unitaria

adf.test(tdd)

    Augmented Dickey-Fuller Test

data:  tdd
Dickey-Fuller = -2.745, Lag order = 7, p-value = 0.2625
alternative hypothesis: stationary

tdd.d<-(diff(tdd))
adf.test(tdd.d)

p-value smaller than printed p-value

    Augmented Dickey-Fuller Test

data:  tdd.d
Dickey-Fuller = -5.9301, Lag order = 7, p-value = 0.01
alternative hypothesis: stationary

7. Usar herramientas ACF,PACF,EACF y/o criterios de Akaike/Bayes para construir propuestas de modelos.

ggAcf(tdd.d)

ggPacf(tdd.d)

eacf(tdd.d)

AR/MA
  0 1 2 3 4 5 6 7 8 9 10 11 12 13
0 x x o x o x o o o o o  o  o  o 
1 o x o x o x o o o o o  o  o  o 
2 x o x o o x o o o o o  o  o  o 
3 x o x x o o o o o o o  o  o  o 
4 o x o o o o o o o o o  o  o  o 
5 o x o x x x o o o o o  o  o  o 
6 x x o o o o o o o o o  o  o  o 
7 x x o o o x o o o o o  o  o  o 

ggtsdisplay(tdd.d)

8. Como parte del diagnóstico del modelo, analizar los residuos y aplicar la prueba Ljung-Box.

fit1<-Arima(tdd.d, order=c(1,1,0), seasonal = c(1,0,0))

fit1

Series: tdd.d 
ARIMA(1,1,0)(1,0,0)[12] 

Coefficients:
          ar1     sar1
      -0.6639  -0.1436
s.e.   0.0388   0.0513

sigma^2 estimated as 0.09411:  log likelihood=-88.17
AIC=182.33   AICc=182.4   BIC=194.11

checkresiduals(fit1)

    Ljung-Box test

data:  Residuals from ARIMA(1,1,0)(1,0,0)[12]
Q* = 120.95, df = 22, p-value = 1.221e-15

Model df: 2.   Total lags used: 24

fit2<-Arima(tdd.d, order=c(2,1,0), seasonal = c(1,0,2))

fit2

Call:
seas(x = tdd, x11 = "")

Coefficients:
       AO1991.Aug         LS1995.Feb         AO1995.Aug  AR-Nonseasonal-01  AR-Nonseasonal-02  
           0.9144             0.8403             0.9726             0.4539             0.2903  
AR-Nonseasonal-03     AR-Seasonal-12  MA-Nonseasonal-01     MA-Seasonal-12  
           0.1097             0.4334             0.8006             0.6729  

checkresiduals(fit2)

fit3<-auto.arima(tdd.d)

fit3

Series: tdd.d 
ARIMA(2,0,3)(0,0,1)[12] with zero mean 

Coefficients:
          ar1     ar2     ma1      ma2     ma3     sma1
      -0.0937  0.6787  -0.232  -0.5449  0.2744  -0.1562
s.e.   0.1307  0.1070   0.133   0.1163  0.0605   0.0498

sigma^2 estimated as 0.0568:  log likelihood=8.46
AIC=-2.92   AICc=-2.62   BIC=24.56

checkresiduals(fit3)

    Ljung-Box test

data:  Residuals from ARIMA(2,0,3)(0,0,1)[12] with zero mean
Q* = 19.407, df = 18, p-value = 0.3672

Model df: 6.   Total lags used: 24

9. Usar la función auto.arima() y compare sus resultados con el inciso anterior.

fit4<-auto.arima(tdd.d, stepwise = FALSE, approximation = FALSE)

fit4

Series: tdd.d 
ARIMA(1,0,2)(1,0,0)[12] with zero mean 

Coefficients:
         ar1      ma1     ma2     sar1
      0.7387  -1.0605  0.3659  -0.1550
s.e.  0.1097   0.1097  0.0497   0.0516

sigma^2 estimated as 0.05681:  log likelihood=7.39
AIC=-4.79   AICc=-4.62   BIC=14.85

checkresiduals(fit4)

    Ljung-Box test

data:  Residuals from ARIMA(1,0,2)(1,0,0)[12] with zero mean
Q* = 24.121, df = 20, p-value = 0.2372

Model df: 4.   Total lags used: 24

10. Presente la ecuación final y describa.

fit5<-Arima(tdd.d, order=c(3,0,5), seasonal = c(0,0,2))
fit5

Series: tdd.d 
ARIMA(3,0,5)(0,0,2)[12] with non-zero mean 

Coefficients:
          ar1     ar2     ar3     ma1      ma2      ma3     ma4     ma5     sma1    sma2
      -1.3015  0.3461  0.6734  1.0451  -0.6105  -0.3802  0.3911  0.0751  -0.1701  0.0702
s.e.   0.1633  0.3220  0.1614  0.1891   0.2714   0.1153  0.0828  0.0784   0.0524  0.0520
         mean
      -0.0035
s.e.   0.0126

sigma^2 estimated as 0.05407:  log likelihood=17.53
AIC=-11.07   AICc=-10.2   BIC=36.06

checkresiduals(fit5)

    Ljung-Box test

data:  Residuals from ARIMA(3,0,5)(0,0,2)[12] with non-zero mean
Q* = 13.886, df = 13, p-value = 0.3819

Model df: 11.   Total lags used: 24

11. Realizar el pronóstico ARIMA para dos años.

getrmse<-function(x,h,...)
{
train.end<-time(x)[length(x)-h]
test.start<-time(x)[length(x)-h+1]
train<-window(x,end=train.end)
test<-window(x,start=test.start)
fit<-Arima(train,...)
fc<-forecast(fit,h=h)
return(accuracy(fc,test)[2,"RMSE"])
}

getrmse(tdd,h=24,order=c(1,0,2),seasonal=c(1,0,0))

[1] 0.4334241

getrmse(tdd,h=24,order=c(2,1,1),seasonal=c(1,1,0))

[1] 0.3770017

getrmse(tdd,h=24,order=c(2,1,0),seasonal=c(1,1,0))

[1] 0.2497599

getrmse(tdd,h=24,order=c(2,1,0),seasonal=c(1,1,2))

[1] 0.3915485

getrmse(tdd,h=24,order=c(2,1,1),seasonal=c(1,1,2))

[1] 0.1926313

getrmse(tdd,h=24,order=c(2,1,1),seasonal=c(1,1,1))

[1] 0.1992385

getrmse(tdd,h=24,order=c(2,1,1),seasonal=c(0,1,1))

[1] 0.2174223

getrmse(tdd,h=24,order=c(3,0,4),seasonal=c(1,1,0))

[1] 0.3284223

getrmse(tdd,h=24,order=c(3,0,4),seasonal=c(1,1,1))

[1] 0.3307197

getrmse(tdd,h=24,order=c(3,0,4),seasonal=c(0,1,2))

[1] 0.3475158

detach("package:bindrcpp", unload=TRUE)

<U+393C><U+3E31>bindrcpp<U+393C><U+3E32> namespace cannot be unloaded:
  namespace <U+393C><U+3E31>bindrcpp<U+393C><U+3E32> is imported by <U+393C><U+3E31>dplyr<U+393C><U+3E32> so cannot be unloaded

getrmse(tdd,h=24,order=c(3,1,4),seasonal=c(1,0,0))

[1] 0.2928143

getrmse(tdd,h=24,order=c(3,0,5),seasonal=c(0,0,2))

[1] 0.1101478

getrmse(tdd,h=24,order=c(3,1,1),seasonal=c(0,0,2))

[1] 0.2924821

fit6<-Arima(tdd, order=c(2,1,1),seasonal=c(1,1,2))
fit6

Series: tdd 
ARIMA(2,1,1)(1,1,2)[12] 

Coefficients:
         ar1     ar2      ma1     sar1     sma1    sma2
      0.5265  0.3396  -0.8141  -0.0228  -1.1079  0.1594
s.e.  0.1018  0.0549   0.0864   0.7714   0.7532  0.7576

sigma^2 estimated as 0.05989:  log likelihood=-16.06
AIC=46.13   AICc=46.44   BIC=73.39

autoplot(forecast(fit6,h=12,type='b',xlab='Time',ylab='Color  Property'))

The non-existent type arguments will be ignored.The non-existent xlab arguments will be ignored.The non-existent ylab arguments will be ignored.

12. Utilizando métodos de pronósticos simples y/o suavizamiento exponencial generar un pronóstico para dos años. Elegir el método que presente la raíz del error medio al cuadrado más pequeño.

Para poder realizar el pronostico de nuestra serie, primeramente se realizo un recorte de los ultimos dos años para despues pronosticar nuestra serie mediante los metodos mean(media), naive y naive estacional.

tdd.c <- window(tdd,start=c(1987,4),end=c(2015,12))
autoplot(tdd.c)

autoplot(tdd.c)+forecast::autolayer(meanf(tdd.c,h=42),PI=FALSE,series="Mean")+forecast::autolayer(rwf(tdd.c,h=42),PI=FALSE,series="Naïve")+forecast::autolayer(rwf(tdd.c,drift=TRUE,h=42),PI=FALSE,series="Drift")+ggtitle("Tasa de Desempleo")+xlab("Day")+ylab("año")+guides(colour=guide_legend(title="Forecast"))

13. Usando el an´alisis de transferencia, modelar observaciones cambios estructurales y/o observaciones atípicas en el contexto ARIMAX (punto extra, opcional)

14. Analizar un VAR, necesita incluir los puntos siguientes

a) Buscar otra variable (incluso más variables) y redactar la teoría que se encuentra detrás de la relación entre variables que propone

En este punto, la variable que usaremos es el IPC, durante mucho tiempo se observó alguna relación entre las tasas de inflación y las tasas de desempleo en algunos países. Así se observó que cuando una nación tenía baja inflación, tendía a experimentar tasas de desempleo más altas, y cuando sus tasas de inflación eran más elevadas, las de desempleo eran menores. Como resultado de esta relación empírica prácticamente se estableció que la sociedad debía elegir entre uno de dos males: inflación o desempleo. ¿Existe algún tipo de correlación entre inflación y empleo? Sí. La más conocida es la llamada curva de Phillips, que sirve para representar la relación inversa entre inflación y desempleo. En teoría, nos indica que a mayor inflación, menor desempleo. Es decir, el mercado laboral se resiente cuando se trata de controlar el IPC.

La teoría se basa en el estudio que Phillips, Solow y Samuelson realizaron en la década de 1960, donde el control de la inflación se tradujo en una contracción económica. La curva de Phillips sirvió (y sigue sirviendo) de excusa para que muchos paíeses mantuviesen unas tasas de inflación elevadas a costa de tener una tasa de paro relativamente baja.

base<-read.csv(file.choose())
plot(base)

plot.ts(base)

b) Evaluar la estacionalidad de la nueva variable.

ts.tdd <- ts(base$tdd, start = c(1987,1), frequency = 12)
ts.IPC <- ts(base$IPC, start = c(1987,1), frequency = 12)

c) Proponga el orden de rezagos óptimo para el VAR y presentar la estimación.

VARselect(base[,1:2],lag.max = (12),type = "const")[["selection"]]

AIC(n)  HQ(n)  SC(n) FPE(n) 
    12      2      2     12 

var1<-VAR(base[,1:2],p=2,type = "const")
summary(var1) ##Verificar que las raíces del polinomio caracteristico sean menor a 1

VAR Estimation Results:
========================= 
Endogenous variables: tdd, IPC 
Deterministic variables: const 
Sample size: 374 
Log Likelihood: -27.96 
Roots of the characteristic polynomial:
    1 0.9788 0.6137 0.2784
Call:
VAR(y = base[, 1:2], p = 2, type = "const")


Estimation results for equation tdd: 
==================================== 
tdd = tdd.l1 + IPC.l1 + tdd.l2 + IPC.l2 + const 

       Estimate Std. Error t value Pr(>|t|)    
tdd.l1  0.69948    0.05002  13.984  < 2e-16 ***
IPC.l1  0.02373    0.03870   0.613    0.540    
tdd.l2  0.27257    0.05001   5.451  9.2e-08 ***
IPC.l2 -0.02351    0.03871  -0.607    0.544    
const   0.08070    0.05052   1.597    0.111    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Residual standard error: 0.244 on 369 degrees of freedom
Multiple R-Squared: 0.9422, Adjusted R-squared: 0.9416 
F-statistic:  1504 on 4 and 369 DF,  p-value: < 2.2e-16 


Estimation results for equation IPC: 
==================================== 
IPC = tdd.l1 + IPC.l1 + tdd.l2 + IPC.l2 + const 

        Estimate Std. Error t value Pr(>|t|)    
tdd.l1  0.009824   0.053711   0.183    0.855    
IPC.l1  1.614819   0.041557  38.858   <2e-16 ***
tdd.l2  0.013188   0.053698   0.246    0.806    
IPC.l2 -0.614899   0.041571 -14.791   <2e-16 ***
const   0.049229   0.054247   0.908    0.365    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Residual standard error: 0.2621 on 369 degrees of freedom
Multiple R-Squared:     1,  Adjusted R-squared:     1 
F-statistic: 2.049e+06 on 4 and 369 DF,  p-value: < 2.2e-16 



Covariance matrix of residuals:
          tdd       IPC
tdd 0.0595567 0.0006551
IPC 0.0006551 0.0686740

Correlation matrix of residuals:
        tdd     IPC
tdd 1.00000 0.01024
IPC 0.01024 1.00000

serial1<-serial.test(var1,lags.pt = (12),type = "PT.asymptotic")
serial1

    Portmanteau Test (asymptotic)

data:  Residuals of VAR object var1
Chi-squared = 165.83, df = 40, p-value < 2.2e-16

var2<-VAR(base[,1:2],p=3,type = "const")
summary(var2)

VAR Estimation Results:
========================= 
Endogenous variables: tdd, IPC 
Deterministic variables: const 
Sample size: 373 
Log Likelihood: -26.093 
Roots of the characteristic polynomial:
    1 0.9789 0.5802 0.3679 0.1162 0.03518
Call:
VAR(y = base[, 1:2], p = 3, type = "const")


Estimation results for equation tdd: 
==================================== 
tdd = tdd.l1 + IPC.l1 + tdd.l2 + IPC.l2 + tdd.l3 + IPC.l3 + const 

       Estimate Std. Error t value Pr(>|t|)    
tdd.l1  0.71276    0.05192  13.729  < 2e-16 ***
IPC.l1 -0.01988    0.04902  -0.405    0.685    
tdd.l2  0.29580    0.06170   4.794 2.38e-06 ***
IPC.l2  0.08757    0.08815   0.994    0.321    
tdd.l3 -0.03713    0.05184  -0.716    0.474    
IPC.l3 -0.06754    0.04904  -1.377    0.169    
const   0.07989    0.05063   1.578    0.115    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Residual standard error: 0.2433 on 366 degrees of freedom
Multiple R-Squared: 0.943,  Adjusted R-squared: 0.942 
F-statistic:  1009 on 6 and 366 DF,  p-value: < 2.2e-16 


Estimation results for equation IPC: 
==================================== 
IPC = tdd.l1 + IPC.l1 + tdd.l2 + IPC.l2 + tdd.l3 + IPC.l3 + const 

        Estimate Std. Error t value Pr(>|t|)    
tdd.l1  0.010877   0.056114   0.194    0.846    
IPC.l1  1.629899   0.052984  30.762  < 2e-16 ***
tdd.l2  0.011784   0.066683   0.177    0.860    
IPC.l2 -0.655301   0.095270  -6.878 2.63e-11 ***
tdd.l3  0.001271   0.056026   0.023    0.982    
IPC.l3  0.025312   0.053001   0.478    0.633    
const   0.050069   0.054720   0.915    0.361    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Residual standard error: 0.263 on 366 degrees of freedom
Multiple R-Squared:     1,  Adjusted R-squared:     1 
F-statistic: 1.347e+06 on 6 and 366 DF,  p-value: < 2.2e-16 



Covariance matrix of residuals:
          tdd       IPC
tdd 0.0592107 0.0006755
IPC 0.0006755 0.0691697

Correlation matrix of residuals:
        tdd     IPC
tdd 1.00000 0.01056
IPC 0.01056 1.00000

var3<-VAR(base[,1:2],p=4,type = "const")
summary(var3)

VAR Estimation Results:
========================= 
Endogenous variables: tdd, IPC 
Deterministic variables: const 
Sample size: 372 
Log Likelihood: -22.284 
Roots of the characteristic polynomial:
    1 0.9744 0.5443 0.5443 0.4683 0.4683 0.4005 0.3985
Call:
VAR(y = base[, 1:2], p = 4, type = "const")


Estimation results for equation tdd: 
==================================== 
tdd = tdd.l1 + IPC.l1 + tdd.l2 + IPC.l2 + tdd.l3 + IPC.l3 + tdd.l4 + IPC.l4 + const 

       Estimate Std. Error t value Pr(>|t|)    
tdd.l1  0.71327    0.05226  13.650  < 2e-16 ***
IPC.l1 -0.01985    0.04902  -0.405   0.6857    
tdd.l2  0.31901    0.06388   4.994 9.22e-07 ***
IPC.l2  0.10536    0.09373   1.124   0.2617    
tdd.l3  0.02524    0.06358   0.397   0.6916    
IPC.l3 -0.11286    0.09385  -1.203   0.2299    
tdd.l4 -0.08824    0.05185  -1.702   0.0897 .  
IPC.l4  0.02755    0.04916   0.560   0.5755    
const   0.08900    0.05086   1.750   0.0809 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Residual standard error: 0.2432 on 363 degrees of freedom
Multiple R-Squared: 0.9434, Adjusted R-squared: 0.9422 
F-statistic:   757 on 8 and 363 DF,  p-value: < 2.2e-16 


Estimation results for equation IPC: 
==================================== 
IPC = tdd.l1 + IPC.l1 + tdd.l2 + IPC.l2 + tdd.l3 + IPC.l3 + tdd.l4 + IPC.l4 + const 

       Estimate Std. Error t value Pr(>|t|)    
tdd.l1  0.02067    0.05629   0.367    0.714    
IPC.l1  1.62756    0.05280  30.826  < 2e-16 ***
tdd.l2 -0.02186    0.06881  -0.318    0.751    
IPC.l2 -0.60169    0.10096  -5.960 5.99e-09 ***
tdd.l3 -0.07389    0.06849  -1.079    0.281    
IPC.l3 -0.10811    0.10109  -1.069    0.286    
tdd.l4  0.10375    0.05585   1.858    0.064 .  
IPC.l4  0.08215    0.05295   1.551    0.122    
const   0.04409    0.05478   0.805    0.421    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Residual standard error: 0.262 on 363 degrees of freedom
Multiple R-Squared:     1,  Adjusted R-squared:     1 
F-statistic: 1.011e+06 on 8 and 363 DF,  p-value: < 2.2e-16 



Covariance matrix of residuals:
        tdd     IPC
tdd 0.05915 0.00109
IPC 0.00109 0.06863

Correlation matrix of residuals:
        tdd     IPC
tdd 1.00000 0.01711
IPC 0.01711 1.00000

serial3<-serial.test(var3,lags.pt = 12,type = "PT.asymptotic")
serial3

    Portmanteau Test (asymptotic)

data:  Residuals of VAR object var3
Chi-squared = 144.06, df = 32, p-value = 3.331e-16

plot(serial3,names="IPC")

plot(serial3,names="tdd")

cusum<-stability(var3,type = "OLS-CUSUM")
plot(cusum)

rec.cusum<-stability(var3,type="Rec-CUSUM")
plot(rec.cusum)

Podemosobservar en las 2 Graficas, que tanto la tasa de desempleo y el IPC se mantienen en una zana contante y que estan dentro de los parametros, tambien observamos que no hay inestabilidad estructural.

base3 <- base2[-c(1:12),]  
plot(base3)

VARselect(base3, lag.max = 12)

$`selection`
AIC(n)  HQ(n)  SC(n) FPE(n) 
    12      1      1     12 

$criteria
                  1            2            3            4            5            6
AIC(n) -5.448043828 -5.435807041 -5.421545652 -5.442583684 -5.437933329 -5.440682462
HQ(n)  -5.421777717 -5.392030188 -5.360258058 -5.363785349 -5.341624253 -5.326862645
SC(n)  -5.382047483 -5.325813131 -5.267554178 -5.244594647 -5.195946728 -5.154698297
FPE(n)  0.004304721  0.004357734  0.004420356  0.004328383  0.004348639  0.004336817
                  7            8            9          10           11           12
AIC(n) -5.443706514 -5.424644215 -5.405233158 -5.42877821 -5.435680639 -5.597025509
HQ(n)  -5.312375956 -5.275802916 -5.238881118 -5.24491543 -5.234307117 -5.378141246
SC(n)  -5.113724786 -5.050664922 -4.987256302 -4.96680379 -4.929708655 -5.047055961
FPE(n)  0.004323879  0.004407302  0.004493959  0.00438971  0.004359907  0.003710672

var12 <- VAR(base3, p=12)
var12

VAR Estimation Results:
======================= 

Estimated coefficients for equation IPC: 
======================================== 
Call:
IPC = IPC.l1 + tdd.l1 + IPC.l2 + tdd.l2 + IPC.l3 + tdd.l3 + IPC.l4 + tdd.l4 + IPC.l5 + tdd.l5 + IPC.l6 + tdd.l6 + IPC.l7 + tdd.l7 + IPC.l8 + tdd.l8 + IPC.l9 + tdd.l9 + IPC.l10 + tdd.l10 + IPC.l11 + tdd.l11 + IPC.l12 + tdd.l12 + const 

      IPC.l1       tdd.l1       IPC.l2       tdd.l2       IPC.l3       tdd.l3       IPC.l4 
 0.409091415  0.060367459 -0.009991214  0.032670331 -0.020740531  0.046580793  0.102588844 
      tdd.l4       IPC.l5       tdd.l5       IPC.l6       tdd.l6       IPC.l7       tdd.l7 
 0.154294127 -0.013610391  0.088414037  0.025181607  0.126693026  0.064188138  0.095266978 
      IPC.l8       tdd.l8       IPC.l9       tdd.l9      IPC.l10      tdd.l10      IPC.l11 
 0.033639028 -0.026548320  0.084550288 -0.107012145 -0.015616196 -0.172109166  0.071161258 
     tdd.l11      IPC.l12      tdd.l12        const 
-0.027032542 -0.411724311  0.005040544  0.006237320 


Estimated coefficients for equation tdd: 
======================================== 
Call:
tdd = IPC.l1 + tdd.l1 + IPC.l2 + tdd.l2 + IPC.l3 + tdd.l3 + IPC.l4 + tdd.l4 + IPC.l5 + tdd.l5 + IPC.l6 + tdd.l6 + IPC.l7 + tdd.l7 + IPC.l8 + tdd.l8 + IPC.l9 + tdd.l9 + IPC.l10 + tdd.l10 + IPC.l11 + tdd.l11 + IPC.l12 + tdd.l12 + const 

       IPC.l1        tdd.l1        IPC.l2        tdd.l2        IPC.l3        tdd.l3 
 0.0490863711 -0.2838683881  0.0718910942 -0.0030042974 -0.0806008764  0.0982968177 
       IPC.l4        tdd.l4        IPC.l5        tdd.l5        IPC.l6        tdd.l6 
 0.1169841804  0.1624157584 -0.1271650999  0.0345887120 -0.0133020398  0.1211837230 
       IPC.l7        tdd.l7        IPC.l8        tdd.l8        IPC.l9        tdd.l9 
 0.0645135903 -0.0542317020 -0.0064755929 -0.0117792907 -0.0442486440  0.0247976959 
      IPC.l10       tdd.l10       IPC.l11       tdd.l11       IPC.l12       tdd.l12 
 0.0369768621  0.0345711688 -0.0243326036 -0.0462242283 -0.0351567562 -0.1225423937 
        const 
-0.0003585096 

d) Probar la correlación serial y normalidad en los residuos.

causality(var12, cause='tdd')$Granger

    Granger causality H0: tdd do not Granger-cause IPC

data:  VAR object var12
F-Test = 2.165, df1 = 12, df2 = 652, p-value = 0.01195

causality(var12, cause='IPC')$Granger

    Granger causality H0: IPC do not Granger-cause tdd

data:  VAR object var12
F-Test = 1.3065, df1 = 12, df2 = 652, p-value = 0.2097

plot(irf(var12, impulse='tdd', response = 'IPC'))

plot(irf(var12, impulse='IPC', response = 'tdd'))

var.serial <- serial.test(var12, lags.pt = 12, type="PT.asymptotic")
var.serial

    Portmanteau Test (asymptotic)

data:  Residuals of VAR object var12
Chi-squared = 33.752, df = 0, p-value < 2.2e-16

var.norm <- normality.test(var12, multivariate.only = T)
var.norm

$`JB`

    JB-Test (multivariate)

data:  Residuals of VAR object var12
Chi-squared = 2072.9, df = 4, p-value < 2.2e-16


$Skewness

    Skewness only (multivariate)

data:  Residuals of VAR object var12
Chi-squared = 145.38, df = 2, p-value < 2.2e-16


$Kurtosis

    Kurtosis only (multivariate)

data:  Residuals of VAR object var12
Chi-squared = 1927.6, df = 2, p-value < 2.2e-16

e) Resultados de la prueba de causalidad de Granger

causality(var3,cause="tdd")$Granger                     

    Granger causality H0: tdd do not Granger-cause IPC

data:  VAR object var3
F-Test = 1.6424, df1 = 4, df2 = 726, p-value = 0.1617

causality(var3,cause="IPC")$Grange

    Granger causality H0: IPC do not Granger-cause tdd

data:  VAR object var3
F-Test = 0.73718, df1 = 4, df2 = 726, p-value = 0.5668

Se puede observar que no se rechaza la hipotesis nula de no Grangercausalidad. ##f) Crear un pronóstico para dos años usando el VAR especificado.

predictions<-predict(var3,n.ahead=10,ci=0.95)
class(predictions)

[1] "varprd"

##  [1]  "varprd"
plot(predictions,names="IPC")

fanchart(predictions,names="tdd")

g) Interpretar las funciones impulso-respuesta.

plot(irf(var3,impulse = "tdd",response = "IPC",ortho=T))

plot(irf(var3,impulse = "IPC",response = "tdd", ortho=T))

15. Evaluar la potencial cointegracióon usando la metodologíia de Engle-Granger y Johansen.

dwtest(lr.reg)

    Durbin-Watson test

data:  lr.reg
DW = 2.5246, p-value = 1
alternative hypothesis: true autocorrelation is greater than 0

ur.df(ts.IPC)

############################################################### 
# Augmented Dickey-Fuller Test Unit Root / Cointegration Test # 
############################################################### 

The value of the test statistic is: 4.8972 

ur.df(ts.tdd)

############################################################### 
# Augmented Dickey-Fuller Test Unit Root / Cointegration Test # 
############################################################### 

The value of the test statistic is: -0.8139 

lr.reg<-dynlm(ts.tdd~L(ts.tdd)+ts.IPC+L(ts.IPC))
summary(lr.reg)

Time series regression with "ts" data:
Start = 1987(2), End = 2018(4)

Call:
dynlm(formula = ts.tdd ~ L(ts.tdd) + ts.IPC + L(ts.IPC))

Residuals:
    Min      1Q  Median      3Q     Max 
-0.7528 -0.1354 -0.0090  0.1471  0.9203 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.11297    0.05207   2.169   0.0307 *  
L(ts.tdd)    0.96313    0.01388  69.412   <2e-16 ***
ts.IPC       0.02213    0.03985   0.555   0.5791    
L(ts.IPC)   -0.02185    0.03986  -0.548   0.5839    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2533 on 371 degrees of freedom
Multiple R-squared:  0.9375,    Adjusted R-squared:  0.937 
F-statistic:  1855 on 3 and 371 DF,  p-value: < 2.2e-16

dwtest(lr.reg)

    Durbin-Watson test

data:  lr.reg
DW = 2.5246, p-value = 1
alternative hypothesis: true autocorrelation is greater than 0

error<-residuals(lr.reg)
plot(error)
abline(h=0)

ur.df(error)

############################################################### 
# Augmented Dickey-Fuller Test Unit Root / Cointegration Test # 
############################################################### 

The value of the test statistic is: -14.4901 

data<-ts.intersect(ts.IPC,ts.tdd,error)
dyn<-dynlm(d(ts.tdd)~L(error)+L(d(ts.IPC))+L(d(ts.tdd)),data=data)
summary(dyn)

Time series regression with "ts" data:
Start = 1987(4), End = 2018(4)

Call:
dynlm(formula = d(ts.tdd) ~ L(error) + L(d(ts.IPC)) + L(d(ts.tdd)), 
    data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.72952 -0.14623 -0.01204  0.12488  0.93426 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
(Intercept)  -0.003095   0.018370  -0.168   0.8663  
L(error)     -0.714139   0.365767  -1.952   0.0516 .
L(d(ts.IPC))  0.005410   0.038084   0.142   0.8871  
L(d(ts.tdd))  0.417836   0.362211   1.154   0.2494  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2431 on 369 degrees of freedom
Multiple R-squared:  0.09021,   Adjusted R-squared:  0.08281 
F-statistic:  12.2 on 3 and 369 DF,  p-value: 1.259e-07