#Bibliotecas
#-----------------------
library(vars)
## Warning: package 'vars' was built under R version 4.3.3
## Loading required package: MASS
## Loading required package: strucchange
## Warning: package 'strucchange' was built under R version 4.3.3
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: sandwich
## Loading required package: urca
## Loading required package: lmtest
library(tseries)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(forecast)
library(urca)
library(highcharter)
## Warning: package 'highcharter' was built under R version 4.3.3
## 
## Attaching package: 'highcharter'
## The following object is masked from 'package:lmtest':
## 
##     unemployment
# Datos
#------------------------------------
data(Canada)
Canada = as.data.frame(Canada)
Canada
##           e     prod       rw     U
## 1  929.6105 405.3665 386.1361  7.53
## 2  929.8040 404.6398 388.1358  7.70
## 3  930.3184 403.8149 390.5401  7.47
## 4  931.4277 404.2158 393.9638  7.27
## 5  932.6620 405.0467 396.7647  7.37
## 6  933.5509 404.4167 400.0217  7.13
## 7  933.5315 402.8191 400.7515  7.40
## 8  933.0769 401.9773 405.7335  8.33
## 9  932.1238 402.0897 409.0504  8.83
## 10 930.6359 401.3067 411.3984 10.43
## 11 929.0971 401.6302 413.0194 12.20
## 12 928.5633 401.5638 415.1670 12.77
## 13 929.0694 402.8157 414.6621 12.43
## 14 930.2655 403.1421 415.7319 12.23
## 15 931.6770 403.0786 416.2315 11.70
## 16 932.1390 403.7188 418.1439 11.20
## 17 932.2767 404.8668 419.7352 11.27
## 18 932.8328 405.6362 420.4842 11.47
## 19 933.7334 405.1363 420.9309 11.30
## 20 934.1772 406.0246 422.1124 11.17
## 21 934.5928 406.4123 423.6278 11.00
## 22 935.6067 406.3009 423.9887 10.63
## 23 936.5111 406.3354 424.1902 10.27
## 24 937.4201 406.7737 426.1270 10.20
## 25 938.4159 405.1525 426.8578  9.67
## 26 938.9992 404.9298 426.7457  9.60
## 27 939.2354 404.5765 426.8858  9.60
## 28 939.6795 404.1995 428.8403  9.50
## 29 940.2497 405.9499 430.1223  9.50
## 30 941.4358 405.8221 430.2307  9.03
## 31 942.2981 406.4463 430.3930  8.70
## 32 943.5322 407.0512 432.0284  8.13
## 33 944.3490 407.9460 433.3886  7.87
## 34 944.8215 408.1796 433.9641  7.67
## 35 945.0671 408.5998 434.4844  7.80
## 36 945.8067 409.0906 436.1569  7.73
## 37 946.8697 408.7042 438.2651  7.57
## 38 946.8766 408.9803 438.7636  7.57
## 39 947.2497 408.3287 439.9498  7.33
## 40 947.6513 407.8857 441.8359  7.57
## 41 948.1840 407.2605 443.1769  7.63
## 42 948.3492 406.7752 444.3592  7.60
## 43 948.0322 406.1794 444.5236  8.17
## 44 947.1065 405.4398 446.9694  9.20
## 45 946.0796 403.2800 450.1586 10.17
## 46 946.1838 403.3649 451.5464 10.33
## 47 946.2258 403.3807 452.2984 10.40
## 48 945.9978 404.0032 453.1201 10.37
## 49 945.5183 404.4774 453.9991 10.60
## 50 945.3514 404.7868 454.9552 11.00
## 51 945.2918 405.2710 455.4824 11.40
## 52 945.4008 405.3830 456.1009 11.73
## 53 945.9058 405.1564 457.2027 11.07
## 54 945.9035 406.4700 457.3886 11.67
## 55 946.3190 406.2293 457.7799 11.47
## 56 946.5796 406.7265 457.5535 11.30
## 57 946.7800 408.5785 458.8024 10.97
## 58 947.6283 409.6767 459.0564 10.63
## 59 948.6221 410.3858 459.1578 10.10
## 60 949.3992 410.5395 459.7037  9.67
## 61 949.9481 410.4453 459.7037  9.53
## 62 949.7945 410.6256 460.0258  9.47
## 63 949.9534 410.8672 461.0257  9.50
## 64 950.2502 411.2359 461.3039  9.27
## 65 950.5380 410.6637 461.4031  9.50
## 66 950.7871 410.8085 462.9277  9.43
## 67 950.8695 412.1160 464.6888  9.70
## 68 950.9281 412.9994 465.0717  9.90
## 69 951.8457 412.9551 464.2851  9.43
## 70 952.6005 412.8241 464.0344  9.30
## 71 953.5976 413.0489 463.4535  8.87
## 72 954.1434 413.6110 465.0717  8.77
## 73 954.5426 413.6048 466.0889  8.60
## 74 955.2631 412.9684 466.6171  8.33
## 75 956.0561 412.2659 465.7478  8.17
## 76 956.7966 412.9106 465.8995  8.03
## 77 957.3865 413.8294 466.4099  7.90
## 78 958.0634 414.2242 466.9552  7.87
## 79 958.7166 415.1678 467.6281  7.53
## 80 959.4881 415.7016 467.7026  6.93
## 81 960.3625 416.8674 469.1348  6.80
## 82 960.7834 417.6104 469.3364  6.70
## 83 961.0290 418.0030 470.0117  6.93
## 84 961.7657 417.2667 469.6472  6.87
# Plot de las series temporales
#--------------------------------
layout(matrix(1:4, nrow = 2, ncol = 2))
plot.ts(Canada$e, main = "Employment", ylab ="", xlab="")
plot.ts(Canada$prod, main = "Productivity", ylab ="", xlab="")
plot.ts(Canada$rw, main = "Real Wage", ylab ="", xlab="")
plot.ts(Canada$U, main = "Unemployment Rate", ylab ="", xlab="")

var_canada <- ts(Canada[,c(1,2,3,4)],frequency = 12)
plot.ts(var_canada)

Estacionariedad

#Empleo
#----------------------------
pp.test(var_canada[,1])
## 
##  Phillips-Perron Unit Root Test
## 
## data:  var_canada[, 1]
## Dickey-Fuller Z(alpha) = -5.8701, Truncation lag parameter = 3, p-value
## = 0.7735
## alternative hypothesis: stationary

P valor = 0.77 > 0.05 –> Acepto H0 –> Raíz Unitaria (No estacionaria)

#Productividad
#-----------------------------
pp.test(var_canada[,2])
## 
##  Phillips-Perron Unit Root Test
## 
## data:  var_canada[, 2]
## Dickey-Fuller Z(alpha) = -7.4247, Truncation lag parameter = 3, p-value
## = 0.6816
## alternative hypothesis: stationary

P valor = 0.68 > 0.05 –> Acepto H0 –> Raíz Unitaria (No estacionaria)

# Salario real
#---------------------
pp.test(var_canada[,3])
## 
##  Phillips-Perron Unit Root Test
## 
## data:  var_canada[, 3]
## Dickey-Fuller Z(alpha) = -4.6296, Truncation lag parameter = 3, p-value
## = 0.8468
## alternative hypothesis: stationary

P valor = 0.84 > 0.05 –> Acepto H0 –> Raíz Unitaria (No estacionaria)

# Tasa de desempleo
#---------------------
pp.test(var_canada[,4])
## 
##  Phillips-Perron Unit Root Test
## 
## data:  var_canada[, 4]
## Dickey-Fuller Z(alpha) = -7.1861, Truncation lag parameter = 3, p-value
## = 0.6957
## alternative hypothesis: stationary

P valor = 0.69 > 0.05 –> Acepto H0 –> Raíz Unitaria (No estacionaria)

# ADF empleo
#-------------------------------
adf.test(var_canada[,1])
## 
##  Augmented Dickey-Fuller Test
## 
## data:  var_canada[, 1]
## Dickey-Fuller = -2.148, Lag order = 4, p-value = 0.5152
## alternative hypothesis: stationary

P valor = 0.51 > 0.05 –> Acepto H0 –> Raíz Unitaria (No estacionaria)

# ADF productividad
#---------------------------------------
adf.test(var_canada[,2])
## 
##  Augmented Dickey-Fuller Test
## 
## data:  var_canada[, 2]
## Dickey-Fuller = -2.8725, Lag order = 4, p-value = 0.218
## alternative hypothesis: stationary

P valor = 0.218 > 0.05 –> Acepto H0 –> Raíz Unitaria (No estacionaria)

# ADF salario real
#---------------------------------------
adf.test(var_canada[,3])
## 
##  Augmented Dickey-Fuller Test
## 
## data:  var_canada[, 3]
## Dickey-Fuller = -2.0558, Lag order = 4, p-value = 0.553
## alternative hypothesis: stationary

P valor = 0.553 > 0.05 –> Acepto H0 –> Raíz Unitaria (No estacionaria)

# ADF tasa de desempleo
#---------------------------------------
adf.test(var_canada[,4])
## 
##  Augmented Dickey-Fuller Test
## 
## data:  var_canada[, 4]
## Dickey-Fuller = -2.5988, Lag order = 4, p-value = 0.3303
## alternative hypothesis: stationary

P valor = 0.3303 > 0.05 –> Acepto H0 –> Raíz Unitaria (No estacionaria)

Las series temporales en niveles no son estaiconarias y debemos estudiar su primera diferencia.

# PP primera diferencia del empleo 
#---------------------------------------
pp.test(diff(var_canada[,1]))
## 
##  Phillips-Perron Unit Root Test
## 
## data:  diff(var_canada[, 1])
## Dickey-Fuller Z(alpha) = -26.227, Truncation lag parameter = 3, p-value
## = 0.01233
## alternative hypothesis: stationary

P valor = 0.01233 > 0.05 –> Rechazo H0 –> NO Raíz Unitaria (Sí estacionaria)

# PP primera diferencia de la productividad 
#---------------------------------------
pp.test(diff(var_canada[,2]))
## Warning in pp.test(diff(var_canada[, 2])): p-value smaller than printed p-value
## 
##  Phillips-Perron Unit Root Test
## 
## data:  diff(var_canada[, 2])
## Dickey-Fuller Z(alpha) = -60.867, Truncation lag parameter = 3, p-value
## = 0.01
## alternative hypothesis: stationary

P valor = 0.01 > 0.05 –> Rechazo H0 –> NO Raíz Unitaria (Sí estacionaria)

# PP primera diferencia del salario real
#---------------------------------------
pp.test(diff(var_canada[,3]))
## Warning in pp.test(diff(var_canada[, 3])): p-value smaller than printed p-value
## 
##  Phillips-Perron Unit Root Test
## 
## data:  diff(var_canada[, 3])
## Dickey-Fuller Z(alpha) = -61.576, Truncation lag parameter = 3, p-value
## = 0.01
## alternative hypothesis: stationary

P valor = 0.01 > 0.05 –> Rechazo H0 –> NO Raíz Unitaria (Sí estacionaria)

# PP primera diferencia de la tasa de desempleo
#---------------------------------------
pp.test(diff(var_canada[,4]))
## Warning in pp.test(diff(var_canada[, 4])): p-value smaller than printed p-value
## 
##  Phillips-Perron Unit Root Test
## 
## data:  diff(var_canada[, 4])
## Dickey-Fuller Z(alpha) = -38.946, Truncation lag parameter = 3, p-value
## = 0.01
## alternative hypothesis: stationary

P valor = 0.01 > 0.05 –> Rechazo H0 –> NO Raíz Unitaria (Sí estacionaria)

# ADF de la primera diferencia del empleo
#--------------------------------------------
adf.test(diff(var_canada[,1]))
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(var_canada[, 1])
## Dickey-Fuller = -3.2668, Lag order = 4, p-value = 0.08274
## alternative hypothesis: stationary

Al 10% sí es estacionaria, pero no al 5%.

# ADF de la primera diferencia de la productividad
#--------------------------------------------
adf.test(diff(var_canada[,2]))
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(var_canada[, 2])
## Dickey-Fuller = -3.2361, Lag order = 4, p-value = 0.08775
## alternative hypothesis: stationary

Al 10% sí es estacionaria, pero no al 5%.

# ADF de la primera diferencia del salario real
#--------------------------------------------
adf.test(diff(var_canada[,3]))
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(var_canada[, 3])
## Dickey-Fuller = -3.2352, Lag order = 4, p-value = 0.08789
## alternative hypothesis: stationary

Al 10% sí es estacionaria, pero no al 5%.

# ADF de la primera diferencia de la tasa de desempleo
#--------------------------------------------
adf.test(diff(var_canada[,4]))
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff(var_canada[, 4])
## Dickey-Fuller = -3.7325, Lag order = 4, p-value = 0.02698
## alternative hypothesis: stationary

Al 10% y al 5% sí es estacionaria.

Con PP: son estacionarias todas las series al 5% Con ADF: Todas las series son estacionarias al 10%.

Diferencia de las variables

# Primera diferencia de todas las variables
#-------------------------------------
var_canada_dif <- diff(var_canada[,c(1,2,3,4)])
head(var_canada_dif)
##                 e       prod        rw     U
## Feb 1  0.19347066 -0.7266325 1.9996500  0.17
## Mar 1  0.51440302 -0.8249509 2.4043538 -0.23
## Apr 1  1.10929985  0.4008901 3.4237043 -0.20
## May 1  1.23431817  0.8309404 2.8008737  0.10
## Jun 1  0.88893413 -0.6299749 3.2570107 -0.24
## Jul 1 -0.01941353 -1.5976113 0.7297971  0.27
#Gráfico de las primeras diferencias
#-----------------------------------------
plot.ts(var_canada_dif)

hchart(var_canada_dif)
## Warning: Deprecated function. Use the `create_axis` function.

Una vez las series son estacionarias, se procede a la estimación

# Selección del orden del VAR
#------------------------------------
VARselect(var_canada_dif, lag.max=8, type = "const")
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      2      1      1      2 
## 
## $criteria
##                   1            2           3            4            5
## AIC(n) -6.270032208 -6.285458074 -6.03818873 -6.035665365 -5.808725091
## HQ(n)  -6.023272860 -5.841291247 -5.39661443 -5.196683581 -4.772335828
## SC(n)  -5.652035378 -5.173063780 -4.43139698 -3.934476142 -3.213138403
## FPE(n)  0.001893667  0.001871884  0.00241986  0.002469804  0.003191486
##                   6            7           8
## AIC(n) -5.558322503 -5.370524857 -5.33131938
## HQ(n)  -4.324525761 -3.939320637 -3.70270768
## SC(n)  -2.468338351 -1.786143241 -1.25254030
## FPE(n)  0.004286004  0.005511821  0.00626064

$criteria 2
AIC(n) -6.285458074 HQ(n) -5.841291247 SC(n) -5.173063780 FPE(n) 0.001871884

Usaremos dos retardos (que son los que minimizan el AIC)

#Modelo VAR p = 2
#----------------------------
modelo_var <- VAR(Canada, p=2)
summary(modelo_var)
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: e, prod, rw, U 
## Deterministic variables: const 
## Sample size: 82 
## Log Likelihood: -175.819 
## Roots of the characteristic polynomial:
## 0.995 0.9081 0.9081 0.7381 0.7381 0.1856 0.1429 0.1429
## Call:
## VAR(y = Canada, p = 2)
## 
## 
## Estimation results for equation e: 
## ================================== 
## e = e.l1 + prod.l1 + rw.l1 + U.l1 + e.l2 + prod.l2 + rw.l2 + U.l2 + const 
## 
##           Estimate Std. Error t value Pr(>|t|)    
## e.l1     1.638e+00  1.500e-01  10.918  < 2e-16 ***
## prod.l1  1.673e-01  6.114e-02   2.736  0.00780 ** 
## rw.l1   -6.312e-02  5.524e-02  -1.143  0.25692    
## U.l1     2.656e-01  2.028e-01   1.310  0.19444    
## e.l2    -4.971e-01  1.595e-01  -3.116  0.00262 ** 
## prod.l2 -1.017e-01  6.607e-02  -1.539  0.12824    
## rw.l2    3.844e-03  5.552e-02   0.069  0.94499    
## U.l2     1.327e-01  2.073e-01   0.640  0.52418    
## const   -1.370e+02  5.585e+01  -2.453  0.01655 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.3628 on 73 degrees of freedom
## Multiple R-Squared: 0.9985,  Adjusted R-squared: 0.9984 
## F-statistic:  6189 on 8 and 73 DF,  p-value: < 2.2e-16 
## 
## 
## Estimation results for equation prod: 
## ===================================== 
## prod = e.l1 + prod.l1 + rw.l1 + U.l1 + e.l2 + prod.l2 + rw.l2 + U.l2 + const 
## 
##           Estimate Std. Error t value Pr(>|t|)    
## e.l1      -0.17277    0.26977  -0.640  0.52390    
## prod.l1    1.15043    0.10995  10.464 3.57e-16 ***
## rw.l1      0.05130    0.09934   0.516  0.60710    
## U.l1      -0.47850    0.36470  -1.312  0.19362    
## e.l2       0.38526    0.28688   1.343  0.18346    
## prod.l2   -0.17241    0.11881  -1.451  0.15104    
## rw.l2     -0.11885    0.09985  -1.190  0.23778    
## U.l2       1.01592    0.37285   2.725  0.00805 ** 
## const   -166.77552  100.43388  -1.661  0.10109    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.6525 on 73 degrees of freedom
## Multiple R-Squared: 0.9787,  Adjusted R-squared: 0.9764 
## F-statistic: 419.3 on 8 and 73 DF,  p-value: < 2.2e-16 
## 
## 
## Estimation results for equation rw: 
## =================================== 
## rw = e.l1 + prod.l1 + rw.l1 + U.l1 + e.l2 + prod.l2 + rw.l2 + U.l2 + const 
## 
##           Estimate Std. Error t value Pr(>|t|)    
## e.l1     -0.268833   0.322619  -0.833    0.407    
## prod.l1  -0.081065   0.131487  -0.617    0.539    
## rw.l1     0.895478   0.118800   7.538 1.04e-10 ***
## U.l1      0.012130   0.436149   0.028    0.978    
## e.l2      0.367849   0.343087   1.072    0.287    
## prod.l2  -0.005181   0.142093  -0.036    0.971    
## rw.l2     0.052677   0.119410   0.441    0.660    
## U.l2     -0.127708   0.445892  -0.286    0.775    
## const   -33.188339 120.110525  -0.276    0.783    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.7803 on 73 degrees of freedom
## Multiple R-Squared: 0.9989,  Adjusted R-squared: 0.9987 
## F-statistic:  8009 on 8 and 73 DF,  p-value: < 2.2e-16 
## 
## 
## Estimation results for equation U: 
## ================================== 
## U = e.l1 + prod.l1 + rw.l1 + U.l1 + e.l2 + prod.l2 + rw.l2 + U.l2 + const 
## 
##          Estimate Std. Error t value Pr(>|t|)    
## e.l1     -0.58076    0.11563  -5.023 3.49e-06 ***
## prod.l1  -0.07812    0.04713  -1.658 0.101682    
## rw.l1     0.01866    0.04258   0.438 0.662463    
## U.l1      0.61893    0.15632   3.959 0.000173 ***
## e.l2      0.40982    0.12296   3.333 0.001352 ** 
## prod.l2   0.05212    0.05093   1.023 0.309513    
## rw.l2     0.04180    0.04280   0.977 0.331928    
## U.l2     -0.07117    0.15981  -0.445 0.657395    
## const   149.78056   43.04810   3.479 0.000851 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.2797 on 73 degrees of freedom
## Multiple R-Squared: 0.9726,  Adjusted R-squared: 0.9696 
## F-statistic:   324 on 8 and 73 DF,  p-value: < 2.2e-16 
## 
## 
## 
## Covariance matrix of residuals:
##              e      prod       rw        U
## e     0.131635 -0.007469 -0.04210 -0.06909
## prod -0.007469  0.425711  0.06461  0.01392
## rw   -0.042099  0.064613  0.60886  0.03422
## U    -0.069087  0.013923  0.03422  0.07821
## 
## Correlation matrix of residuals:
##             e     prod      rw       U
## e     1.00000 -0.03155 -0.1487 -0.6809
## prod -0.03155  1.00000  0.1269  0.0763
## rw   -0.14870  0.12691  1.0000  0.1568
## U    -0.68090  0.07630  0.1568  1.0000

Diagnosis

roots(modelo_var)
## [1] 0.9950338 0.9081062 0.9081062 0.7380565 0.7380565 0.1856381 0.1428889
## [8] 0.1428889

Al ser inferiores a uno todas las raíces del polinomio, el VAR es estacionario.

Contraste de autocorrelación de los errores

serial.test(modelo_var,lags.pt=16, type ="PT.adjusted")
## 
##  Portmanteau Test (adjusted)
## 
## data:  Residuals of VAR object modelo_var
## Chi-squared = 231.59, df = 224, p-value = 0.3497

El p valor es 0.3497 > 0.05. Aceptamos H0, es decir, que los errores no están correlacionados.

Prueba de Normaliad

normality.test(modelo_var, multivariate.only = FALSE)
## $e
## 
##  JB-Test (univariate)
## 
## data:  Residual of e equation
## Chi-squared = 0.15347, df = 2, p-value = 0.9261
## 
## 
## $prod
## 
##  JB-Test (univariate)
## 
## data:  Residual of prod equation
## Chi-squared = 4.2651, df = 2, p-value = 0.1185
## 
## 
## $rw
## 
##  JB-Test (univariate)
## 
## data:  Residual of rw equation
## Chi-squared = 0.33482, df = 2, p-value = 0.8459
## 
## 
## $U
## 
##  JB-Test (univariate)
## 
## data:  Residual of U equation
## Chi-squared = 0.56643, df = 2, p-value = 0.7534
## 
## 
## $JB
## 
##  JB-Test (multivariate)
## 
## data:  Residuals of VAR object modelo_var
## Chi-squared = 5.094, df = 8, p-value = 0.7475
## 
## 
## $Skewness
## 
##  Skewness only (multivariate)
## 
## data:  Residuals of VAR object modelo_var
## Chi-squared = 1.7761, df = 4, p-value = 0.7769
## 
## 
## $Kurtosis
## 
##  Kurtosis only (multivariate)
## 
## data:  Residuals of VAR object modelo_var
## Chi-squared = 3.3179, df = 4, p-value = 0.5061

Todos los errores son normales, al tener un p valor mayor a un 5% en todos los casos.