library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.1.3
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.5 v dplyr 1.0.7
## v tidyr 1.1.4 v stringr 1.4.0
## v readr 2.1.2 v forcats 0.5.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(lubridate)
##
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union
library(car)
## Warning: package 'car' was built under R version 4.1.3
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
##
## recode
## The following object is masked from 'package:purrr':
##
## some
library(urca)
library(tseries)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
library(aTSA)
##
## Attaching package: 'aTSA'
## The following objects are masked from 'package:tseries':
##
## adf.test, kpss.test, pp.test
## The following object is masked from 'package:graphics':
##
## identify
library(forecast)
## Warning: package 'forecast' was built under R version 4.1.3
##
## Attaching package: 'forecast'
## The following object is masked from 'package:aTSA':
##
## forecast
library(foreign)
library(timsac)
library(vars)
## Warning: package 'vars' was built under R version 4.1.3
## Loading required package: MASS
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
## Loading required package: strucchange
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: sandwich
##
## Attaching package: 'strucchange'
## The following object is masked from 'package:stringr':
##
## boundary
## Loading required package: lmtest
##
## Attaching package: 'vars'
## The following object is masked from 'package:aTSA':
##
## arch.test
library(mFilter)
## Warning: package 'mFilter' was built under R version 4.1.3
library(dynlm)
## Warning: package 'dynlm' was built under R version 4.1.3
library(readxl)
## Warning: package 'readxl' was built under R version 4.1.3
coint<-read_excel(file.choose())
Cointegración_en_R <- read_excel(file.choose())
View(Cointegración_en_R)
coint=as.data.frame(coint)
attach(coint)
class(coint)
## [1] "data.frame"
plot(coint)
plot(PCE)
names(coint)
## [1] "Year" "Quarter" "tiempo" "DPI" "GDP" "PCE" "CP"
## [8] "DIVIDEND"
class(coint)
## [1] "data.frame"
#Generar logaritmos
lnPCE = log(PCE)
lnDPI = log(DPI)
#Crear variables de series de tiempo
DPI.ts = ts(lnDPI, start=c(1974,1), end=c(2007,4), frequency = 4)
PCE.ts = ts(lnPCE, start=c(1974,1), end=c(2007,4), frequency = 4)
ts.plot(DPI.ts)
ts.plot(PCE.ts)
datos1=cbind(DPI.ts, PCE.ts)
plot(cbind(DPI.ts, PCE.ts))
##Cointegracion
#Generar modelo
cor(PCE,DPI)
## [1] 0.9984923
datos1 <- cbind(PCE,DPI)
colnames(datos1) <- cbind( "PCE", "DPI")
lagselect <- VARselect(datos1, lag.max = 8, type = "both")
lagselect$selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 4 2 2 4
lagselect$criteria
## 1 2 3 4 5
## AIC(n) 13.45290 13.31812 13.32334 13.30588 13.32914
## HQ(n) 13.50023 13.38911 13.41800 13.42421 13.47114
## SC(n) 13.57032 13.49424 13.55818 13.59942 13.68140
## FPE(n) 695859.79225 608126.06456 611329.67769 600776.74404 614962.65296
## 6 7 8
## AIC(n) 13.34504 13.36627 13.39234
## HQ(n) 13.51070 13.55560 13.60533
## SC(n) 13.75600 13.83594 13.92072
## FPE(n) 624879.96422 638379.79493 655353.77622
modelo1 <- VAR(datos1, p = 2, season = NULL, exog = NULL, type = "const")
modelo1 = lm(PCE.ts ~ DPI.ts)
summary(modelo1)
##
## Call:
## lm(formula = PCE.ts ~ DPI.ts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.031067 -0.009396 -0.001756 0.007100 0.046117
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.194408 0.021496 9.044 1.51e-15 ***
## DPI.ts 0.960252 0.002806 342.210 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0128 on 134 degrees of freedom
## Multiple R-squared: 0.9989, Adjusted R-squared: 0.9988
## F-statistic: 1.171e+05 on 1 and 134 DF, p-value: < 2.2e-16
residuales=modelo1$residuals
summary(residuales)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.031067 -0.009396 -0.001756 0.000000 0.007100 0.046117
residualPlot(modelo1)
#prueba de raiz unitaria
adf.test(residuales)
## Augmented Dickey-Fuller Test
## alternative: stationary
##
## Type 1: no drift no trend
## lag ADF p.value
## [1,] 0 -6.05 0.01
## [2,] 1 -5.01 0.01
## [3,] 2 -4.90 0.01
## [4,] 3 -5.36 0.01
## [5,] 4 -4.10 0.01
## Type 2: with drift no trend
## lag ADF p.value
## [1,] 0 -6.02 0.01
## [2,] 1 -5.00 0.01
## [3,] 2 -4.89 0.01
## [4,] 3 -5.36 0.01
## [5,] 4 -4.12 0.01
## Type 3: with drift and trend
## lag ADF p.value
## [1,] 0 -6.00 0.01
## [2,] 1 -5.02 0.01
## [3,] 2 -4.93 0.01
## [4,] 3 -5.45 0.01
## [5,] 4 -4.26 0.01
## ----
## Note: in fact, p.value = 0.01 means p.value <= 0.01
y = ur.df(residuales, type = "trend", selectlags = "AIC")
summary(y)
##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression trend
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.038950 -0.005109 0.000418 0.004343 0.039482
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.662e-03 1.724e-03 -0.964 0.3369
## z.lag.1 -3.753e-01 7.480e-02 -5.017 1.69e-06 ***
## tt 1.829e-05 2.192e-05 0.835 0.4055
## z.diff.lag -1.680e-01 8.090e-02 -2.077 0.0398 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.009808 on 130 degrees of freedom
## Multiple R-squared: 0.2775, Adjusted R-squared: 0.2608
## F-statistic: 16.64 on 3 and 130 DF, p-value: 3.279e-09
##
##
## Value of test-statistic is: -5.017 8.6193 12.8112
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -3.99 -3.43 -3.13
## phi2 6.22 4.75 4.07
## phi3 8.43 6.49 5.47
y@teststat
## tau3 phi2 phi3
## statistic -5.016975 8.619292 12.81124
y@cval
## 1pct 5pct 10pct
## tau3 -3.99 -3.43 -3.13
## phi2 6.22 4.75 4.07
## phi3 8.43 6.49 5.47
y2 = ur.df(residuales, type = "drift", selectlags = "AIC")
summary(y2)
##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression drift
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.040038 -0.004852 0.000548 0.004892 0.038497
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0004084 0.0008463 -0.483 0.6302
## z.lag.1 -0.3732403 0.0746724 -4.998 1.82e-06 ***
## z.diff.lag -0.1683817 0.0808055 -2.084 0.0391 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.009797 on 131 degrees of freedom
## Multiple R-squared: 0.2736, Adjusted R-squared: 0.2625
## F-statistic: 24.67 on 2 and 131 DF, p-value: 8.05e-10
##
##
## Value of test-statistic is: -4.9984 12.6098
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1 6.52 4.63 3.81
#si el valor de t es menor en valor absoluto al valor critico del
#los residuos son No estacionarios
y3 = ur.df(residuales, type = "none", selectlags = "AIC")
summary(y3)
##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression none
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.040447 -0.005259 0.000139 0.004482 0.038084
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 -0.37326 0.07446 -5.013 1.69e-06 ***
## z.diff.lag -0.16819 0.08057 -2.088 0.0388 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.009768 on 132 degrees of freedom
## Multiple R-squared: 0.2732, Adjusted R-squared: 0.2622
## F-statistic: 24.81 on 2 and 132 DF, p-value: 7.139e-10
##
##
## Value of test-statistic is: -5.0133
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.58 -1.95 -1.62
#Philips y Oularis
prueba.P0 = ca.po(datos1, type="Pz")
summary(prueba.P0)
##
## ########################################
## # Phillips and Ouliaris Unit Root Test #
## ########################################
##
## Test of type Pz
## detrending of series none
##
## Response PCE :
##
## Call:
## lm(formula = PCE ~ zr - 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -105.869 -11.242 0.804 13.818 60.810
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## zrPCE 0.995909 0.010760 92.559 <2e-16 ***
## zrDPI 0.011342 0.009885 1.147 0.252
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 22.5 on 241 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 4.157e+06 on 2 and 241 DF, p-value: < 2.2e-16
##
##
## Response DPI :
##
## Call:
## lm(formula = DPI ~ zr - 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -163.399 -13.419 4.837 19.280 143.173
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## zrPCE 0.02016 0.01911 1.055 0.293
## zrDPI 0.98901 0.01756 56.320 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 39.96 on 241 degrees of freedom
## Multiple R-squared: 0.9999, Adjusted R-squared: 0.9999
## F-statistic: 1.558e+06 on 2 and 241 DF, p-value: < 2.2e-16
##
##
##
## Value of test-statistic is: 8.912
##
## Critical values of Pz are:
## 10pct 5pct 1pct
## critical values 33.9267 40.8217 55.1911
prueba.P2 = ca.po(datos1, type="Pu")
summary(prueba.P2)
##
## ########################################
## # Phillips and Ouliaris Unit Root Test #
## ########################################
##
## Test of type Pu
## detrending of series none
##
##
## Call:
## lm(formula = z[, 1] ~ z[, -1] - 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -292.36 -124.00 -42.27 2.22 395.85
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z[, -1] 0.918864 0.001929 476.3 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 136.1 on 243 degrees of freedom
## Multiple R-squared: 0.9989, Adjusted R-squared: 0.9989
## F-statistic: 2.269e+05 on 1 and 243 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: 6.3268
##
## Critical values of Pu are:
## 10pct 5pct 1pct
## critical values 20.3933 25.9711 38.3413
#Segun la prueba de Philips y Oularis el modelo presenta datos estacionales
##Corrección de errores
#Generación de tendencia si no están cointegrados
tendencia=seq_along(PCE.ts)
print(tendencia)
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
## [19] 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
## [37] 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54
## [55] 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72
## [73] 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90
## [91] 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108
## [109] 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126
## [127] 127 128 129 130 131 132 133 134 135 136
modelo2=lm(PCE.ts~tendencia+DPI.ts)
residuales2 = modelo2$residuals
residualPlot(modelo2)
DFresiduales2 = ur.df(residuales2, type = "trend", selectlags = "AIC")
summary(DFresiduales2)
##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression trend
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.033155 -0.005175 -0.000228 0.004714 0.040952
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.428e-03 1.616e-03 -0.884 0.3783
## z.lag.1 -3.483e-01 7.237e-02 -4.813 4.06e-06 ***
## tt 1.423e-05 2.054e-05 0.693 0.4898
## z.diff.lag -1.844e-01 8.088e-02 -2.279 0.0243 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.009195 on 130 degrees of freedom
## Multiple R-squared: 0.2688, Adjusted R-squared: 0.2519
## F-statistic: 15.93 on 3 and 130 DF, p-value: 7.031e-09
##
##
## Value of test-statistic is: -4.8128 7.9497 11.7426
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -3.99 -3.43 -3.13
## phi2 6.22 4.75 4.07
## phi3 8.43 6.49 5.47
dlnPCE = diff(PCE.ts)
dlnDPI = diff(DPI.ts)
modelo3 = lm(dlnPCE~dlnDPI)
summary(modelo3)
##
## Call:
## lm(formula = dlnPCE ~ dlnDPI)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.039263 -0.004656 0.000486 0.004847 0.039199
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.005772 0.001040 5.552 1.48e-07 ***
## dlnDPI 0.336966 0.069578 4.843 3.50e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.009362 on 133 degrees of freedom
## Multiple R-squared: 0.1499, Adjusted R-squared: 0.1435
## F-statistic: 23.45 on 1 and 133 DF, p-value: 3.499e-06
res3=modelo3$residuals
res3_1 = lag(res3)
MCE = lm(dlnPCE ~ dlnDPI + res3_1)
summary(MCE)
##
## Call:
## lm(formula = dlnPCE ~ dlnDPI + res3_1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.035477 -0.005194 0.000119 0.005468 0.041716
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.004802 0.001045 4.593 1.01e-05 ***
## dlnDPI 0.423487 0.071830 5.896 2.99e-08 ***
## res3_1 -0.238421 0.087117 -2.737 0.00707 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.009041 on 131 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.216, Adjusted R-squared: 0.204
## F-statistic: 18.05 on 2 and 131 DF, p-value: 1.195e-07
#res debe ser negativo y valor absoluto menor a 1