Taller 3
LUIS SUAREZ
2025-06-21
### Cargamos librerias necesarias
library(readxl)
library(car)## Warning: package 'car' was built under R version 4.3.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.3.3library(tseries)## Warning: package 'tseries' was built under R version 4.3.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(forecast)## Warning: package 'forecast' was built under R version 4.3.3library(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## Warning: package 'zoo' was built under R version 4.3.3##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: sandwich## Warning: package 'sandwich' was built under R version 4.3.3## Loading required package: urca## Warning: package 'urca' was built under R version 4.3.3## Loading required package: lmtest## Warning: package 'lmtest' was built under R version 4.3.3library(urca)
library(tseries)
library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:MASS':
##
## select## The following object is masked from 'package:car':
##
## recode## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union### Carga de datos
datos <- read_excel("DATOS_TALLER_3.xlsx")# Transformamos los valores
datos <- datos %>%
mutate(
Vidrio_manufactura = Vidrio_manufactura / 1000,
Aluminio_manufacturas = Aluminio_manufacturas / 1000
)
attach(datos)
head(datos)## # A tibble: 6 × 5
## MES Vidrio_manufactura Aluminio_manufacturas SECUESTRO EXTORSION
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Enero 16.6 29.3 20 652
## 2 Febrero 22.3 39.2 15 650
## 3 Marzo 21.0 34.5 18 574
## 4 Abril 20.8 34.9 9 595
## 5 Mayo 19.5 34.4 14 609
## 6 Junio 20.8 31.4 17 549### aplicamos logaritmos a las variables
logVM<-log(Vidrio_manufactura)
logAM<-log(Aluminio_manufacturas)
logSEC<-log(SECUESTRO)
logEXT<-log(EXTORSION)#Se crean las series de tiempo
logVM_ts <- ts(logVM, start = c(2015, 1), frequency = 12)
logAM_ts <- ts(logAM, start = c(2015, 1), frequency = 12)
logSEC_ts <- ts(logSEC, start = c(2015, 1), frequency = 12)
logEXT_ts <- ts(logEXT, start = c(2015, 1), frequency = 12)### Reviso la tendencia
ts.plot(cbind(logVM_ts , logAM_ts, logSEC_ts, logEXT_ts),
col = c("blue", "darkgreen", "red", "orange"),
lty = 1, lwd = 2,
main = "Series en logaritmos",
ylab = "Log(valor)")
legend("topleft",
legend = c("logVM", "logAM", "logSEC", "logEXT"),
col = c("blue", "darkgreen", "red", "orange"),
lty = 1, bty = "n")
#Verificar la estacionariedad individual de mis variables
# Vidrio_manufactura
adf1_vm <- summary(ur.df(logVM_ts, lags = 1))
adf1_vm #No es estacionaria##
## ###############################################
## # 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.62509 -0.06792 0.00069 0.07449 0.44687
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 9.717e-05 4.280e-03 0.023 0.98192
## z.diff.lag -2.823e-01 8.755e-02 -3.225 0.00164 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.14 on 116 degrees of freedom
## Multiple R-squared: 0.08242, Adjusted R-squared: 0.0666
## F-statistic: 5.21 on 2 and 116 DF, p-value: 0.006813
##
##
## Value of test-statistic is: 0.0227
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.58 -1.95 -1.62adf2_vm <- summary(ur.df(logVM_ts, type = "drift", lags = 1))
adf2_vm #No es estacionaria ##
## ###############################################
## # 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.65398 -0.06960 -0.00492 0.07122 0.36762
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.25471 0.13226 1.926 0.05659 .
## z.lag.1 -0.08407 0.04391 -1.915 0.05802 .
## z.diff.lag -0.23991 0.08930 -2.686 0.00829 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1384 on 115 degrees of freedom
## Multiple R-squared: 0.111, Adjusted R-squared: 0.09553
## F-statistic: 7.179 on 2 and 115 DF, p-value: 0.001154
##
##
## Value of test-statistic is: -1.9147 1.8547
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1 6.52 4.63 3.81adf3_vm<- summary(ur.df(logVM_ts, type = "trend", lags = 1))
adf3_vm #No estacionaria ##
## ###############################################
## # 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.66977 -0.06191 0.00684 0.06641 0.32382
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3407539 0.1406757 2.422 0.0170 *
## z.lag.1 -0.1274814 0.0505286 -2.523 0.0130 *
## tt 0.0007294 0.0004304 1.695 0.0929 .
## z.diff.lag -0.2194572 0.0894012 -2.455 0.0156 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1373 on 114 degrees of freedom
## Multiple R-squared: 0.1328, Adjusted R-squared: 0.11
## F-statistic: 5.821 on 3 and 114 DF, p-value: 0.0009756
##
##
## Value of test-statistic is: -2.523 2.214 3.2989
##
## 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# Aluminio_manufacturas
adf1_am <- summary(ur.df(logAM_ts , lags = 1))
adf1_am #No es estacionaria##
## ###############################################
## # 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.45049 -0.10145 -0.01124 0.08462 0.43567
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 0.0005381 0.0039057 0.138 0.891
## z.diff.lag -0.4077816 0.0836989 -4.872 3.53e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.154 on 116 degrees of freedom
## Multiple R-squared: 0.17, Adjusted R-squared: 0.1556
## F-statistic: 11.88 on 2 and 116 DF, p-value: 2.032e-05
##
##
## Value of test-statistic is: 0.1378
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.58 -1.95 -1.62adf2_am <- summary(ur.df(logAM_ts , type = "drift", lags = 1))
adf2_am # No estacionaria##
## ###############################################
## # 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.48433 -0.10762 -0.01538 0.08675 0.46439
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.49880 0.21202 2.353 0.020339 *
## z.lag.1 -0.13650 0.05837 -2.338 0.021095 *
## z.diff.lag -0.33840 0.08724 -3.879 0.000175 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1511 on 115 degrees of freedom
## Multiple R-squared: 0.208, Adjusted R-squared: 0.1942
## F-statistic: 15.1 on 2 and 115 DF, p-value: 1.507e-06
##
##
## Value of test-statistic is: -2.3384 2.7774
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1 6.52 4.63 3.81adf3_am <- summary(ur.df(logAM_ts , type = "trend", lags = 1))
adf3_am #No Estacionaria ##
## ###############################################
## # 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.50704 -0.09301 -0.01218 0.08751 0.45405
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.7488088 0.2391666 3.131 0.00221 **
## z.lag.1 -0.2230572 0.0702556 -3.175 0.00193 **
## tt 0.0010533 0.0004915 2.143 0.03425 *
## z.diff.lag -0.2962164 0.0881397 -3.361 0.00106 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1488 on 114 degrees of freedom
## Multiple R-squared: 0.2386, Adjusted R-squared: 0.2186
## F-statistic: 11.91 on 3 and 114 DF, p-value: 7.669e-07
##
##
## Value of test-statistic is: -3.1749 3.4401 5.1154
##
## 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# SECUESTRO
adf1_sec <- summary(ur.df(logSEC_ts , lags = 1))
adf1_sec #No es estacionaria##
## ###############################################
## # 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
## -1.89944 -0.24807 0.07132 0.32047 1.41516
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 -0.006028 0.015347 -0.393 0.695
## z.diff.lag -0.444259 0.083579 -5.315 5.22e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4384 on 116 degrees of freedom
## Multiple R-squared: 0.2011, Adjusted R-squared: 0.1873
## F-statistic: 14.6 on 2 and 116 DF, p-value: 2.219e-06
##
##
## Value of test-statistic is: -0.3928
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.58 -1.95 -1.62adf2_sec <- summary(ur.df(logSEC_ts , type = "drift", lags = 1))
adf2_sec #Es estacionaria al 1% 5% y 10%##
## ###############################################
## # 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
## -1.92622 -0.20452 0.01719 0.28550 0.83209
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.08036 0.25701 4.204 5.22e-05 ***
## z.lag.1 -0.41240 0.09773 -4.220 4.91e-05 ***
## z.diff.lag -0.23416 0.09277 -2.524 0.013 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4099 on 115 degrees of freedom
## Multiple R-squared: 0.3074, Adjusted R-squared: 0.2953
## F-statistic: 25.52 on 2 and 115 DF, p-value: 6.734e-10
##
##
## Value of test-statistic is: -4.2196 8.9229
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1 6.52 4.63 3.81adf3_sec <- summary(ur.df(logSEC_ts , type = "trend", lags = 1))
adf3_sec #Es estacionaria al 1% 5% y 10%##
## ###############################################
## # 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
## -1.93041 -0.22209 0.01455 0.28165 0.76139
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.084890 0.254538 4.262 4.19e-05 ***
## z.lag.1 -0.462067 0.100621 -4.592 1.14e-05 ***
## tt 0.002061 0.001142 1.806 0.0736 .
## z.diff.lag -0.214880 0.092488 -2.323 0.0219 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.406 on 114 degrees of freedom
## Multiple R-squared: 0.3266, Adjusted R-squared: 0.3089
## F-statistic: 18.43 on 3 and 114 DF, p-value: 8.089e-10
##
##
## Value of test-statistic is: -4.5922 7.1521 10.7075
##
## 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# EXTORSION
adf1_ext <- summary(ur.df(logEXT_ts , lags = 1))
adf1_ext #No es estacionaria ##
## ###############################################
## # 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.84626 -0.10710 0.01548 0.10592 1.82748
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 0.0002137 0.0041312 0.052 0.959
## z.diff.lag -0.1520081 0.0918887 -1.654 0.101
##
## Residual standard error: 0.2894 on 116 degrees of freedom
## Multiple R-squared: 0.02305, Adjusted R-squared: 0.006208
## F-statistic: 1.369 on 2 and 116 DF, p-value: 0.2586
##
##
## Value of test-statistic is: 0.0517
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.58 -1.95 -1.62adf2_ext <- summary(ur.df(logEXT_ts , type = "drift", lags = 1))
adf2_ext #Es estacionaria 1% 5% y 10%##
## ###############################################
## # 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
## -1.05301 -0.09054 0.01701 0.11895 1.53445
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.22212 0.39505 3.094 0.00248 **
## z.lag.1 -0.18887 0.06125 -3.084 0.00256 **
## z.diff.lag -0.05215 0.09436 -0.553 0.58158
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2792 on 115 degrees of freedom
## Multiple R-squared: 0.09777, Adjusted R-squared: 0.08208
## F-statistic: 6.231 on 2 and 115 DF, p-value: 0.002696
##
##
## Value of test-statistic is: -3.0835 4.7865
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1 6.52 4.63 3.81adf3_ext <- summary(ur.df(logEXT_ts , type = "trend", lags = 1))
adf3_ext #Es estacionaria 1% 5% y 10%##
## ###############################################
## # 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
## -1.13665 -0.09117 0.01782 0.10652 1.28830
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.904771 0.508801 5.709 9.16e-08 ***
## z.lag.1 -0.497490 0.086280 -5.766 7.06e-08 ***
## tt 0.005020 0.001064 4.719 6.80e-06 ***
## z.diff.lag 0.087038 0.091570 0.951 0.344
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2565 on 114 degrees of freedom
## Multiple R-squared: 0.2452, Adjusted R-squared: 0.2253
## F-statistic: 12.34 on 3 and 114 DF, p-value: 4.742e-07
##
##
## Value of test-statistic is: -5.766 11.2026 16.7654
##
## 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### Estimar el modelo VAR
#Hallar diferencias
library(forecast)
ndiffs(logVM_ts) #1## [1] 1ndiffs(logAM_ts) #1## [1] 1ndiffs(logSEC_ts ) #1## [1] 1ndiffs(logEXT_ts) #1## [1] 1#Aplico las diferencias
dlogVM<-diff(logVM_ts,1)
dlogAM <-diff(logAM_ts ,1)
dlogSEC<-diff(logSEC_ts,1)
dlogEXT<-diff(logEXT_ts,1)###Creamos un data frame
datos1 <- cbind(dlogVM,dlogAM,dlogSEC,dlogEXT)
head(datos1)## dlogVM dlogAM dlogSEC dlogEXT
## Feb 2015 0.293667178 0.29215219 -0.28768207 -0.003072199
## Mar 2015 -0.060467302 -0.12743975 0.18232156 -0.124342967
## Apr 2015 -0.009721001 0.01266799 -0.69314718 0.035932009
## May 2015 -0.064884307 -0.01686631 0.44183275 0.023256862
## Jun 2015 0.067125131 -0.09134746 0.19415601 -0.103719826
## Jul 2015 0.042779312 0.12398503 0.05715841 -0.018382871dim(datos1)## [1] 119 4VARselect(datos1, lag.max = 12) ## $selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 3 1 1 3
##
## $criteria
## 1 2 3 4 5
## AIC(n) -1.226224e+01 -1.236025e+01 -1.237893e+01 -1.226240e+01 -1.215604e+01
## HQ(n) -1.205971e+01 -1.199570e+01 -1.185236e+01 -1.157380e+01 -1.130542e+01
## SC(n) -1.176264e+01 -1.146098e+01 -1.107999e+01 -1.056378e+01 -1.005775e+01
## FPE(n) 4.728210e-06 4.292434e-06 4.226614e-06 4.777774e-06 5.366136e-06
## 6 7 8 9 10
## AIC(n) -1.211311e+01 -1.203036e+01 -1.206979e+01 -1.204434e+01 -1.187479e+01
## HQ(n) -1.110046e+01 -1.085570e+01 -1.073310e+01 -1.054562e+01 -1.021406e+01
## SC(n) -9.615135e+00 -9.132716e+00 -8.772466e+00 -8.347340e+00 -7.778120e+00
## FPE(n) 5.683554e-06 6.300998e-06 6.226052e-06 6.619967e-06 8.211438e-06
## 11 12
## AIC(n) -1.192208e+01 -1.184337e+01
## HQ(n) -1.009932e+01 -9.858585e+00
## SC(n) -7.425732e+00 -6.947345e+00
## FPE(n) 8.297297e-06 9.643104e-06###
Var1<-VAR(datos1, p=1)
Var1##
## VAR Estimation Results:
## =======================
##
## Estimated coefficients for equation dlogVM:
## ===========================================
## Call:
## dlogVM = dlogVM.l1 + dlogAM.l1 + dlogSEC.l1 + dlogEXT.l1 + const
##
## dlogVM.l1 dlogAM.l1 dlogSEC.l1 dlogEXT.l1 const
## -0.276804369 0.078506492 0.056140321 -0.019998754 0.002027705
##
##
## Estimated coefficients for equation dlogAM:
## ===========================================
## Call:
## dlogAM = dlogVM.l1 + dlogAM.l1 + dlogSEC.l1 + dlogEXT.l1 + const
##
## dlogVM.l1 dlogAM.l1 dlogSEC.l1 dlogEXT.l1 const
## 0.149950932 -0.448814468 0.026738505 0.033858216 0.003339679
##
##
## Estimated coefficients for equation dlogSEC:
## ============================================
## Call:
## dlogSEC = dlogVM.l1 + dlogAM.l1 + dlogSEC.l1 + dlogEXT.l1 + const
##
## dlogVM.l1 dlogAM.l1 dlogSEC.l1 dlogEXT.l1 const
## -0.069045425 -0.111399931 -0.450897806 0.064890992 0.008070302
##
##
## Estimated coefficients for equation dlogEXT:
## ============================================
## Call:
## dlogEXT = dlogVM.l1 + dlogAM.l1 + dlogSEC.l1 + dlogEXT.l1 + const
##
## dlogVM.l1 dlogAM.l1 dlogSEC.l1 dlogEXT.l1 const
## 0.063536685 -0.179834722 0.024310610 -0.168274124 0.007167226summary(Var1)##
## VAR Estimation Results:
## =========================
## Endogenous variables: dlogVM, dlogAM, dlogSEC, dlogEXT
## Deterministic variables: const
## Sample size: 118
## Log Likelihood: 50.723
## Roots of the characteristic polynomial:
## 0.4876 0.4105 0.2542 0.1925
## Call:
## VAR(y = datos1, p = 1)
##
##
## Estimation results for equation dlogVM:
## =======================================
## dlogVM = dlogVM.l1 + dlogAM.l1 + dlogSEC.l1 + dlogEXT.l1 + const
##
## Estimate Std. Error t value Pr(>|t|)
## dlogVM.l1 -0.276804 0.092652 -2.988 0.00345 **
## dlogAM.l1 0.078506 0.079769 0.984 0.32713
## dlogSEC.l1 0.056140 0.027129 2.069 0.04079 *
## dlogEXT.l1 -0.019999 0.044918 -0.445 0.65701
## const 0.002028 0.012718 0.159 0.87361
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.138 on 113 degrees of freedom
## Multiple R-Squared: 0.1309, Adjusted R-squared: 0.1001
## F-statistic: 4.254 on 4 and 113 DF, p-value: 0.003033
##
##
## Estimation results for equation dlogAM:
## =======================================
## dlogAM = dlogVM.l1 + dlogAM.l1 + dlogSEC.l1 + dlogEXT.l1 + const
##
## Estimate Std. Error t value Pr(>|t|)
## dlogVM.l1 0.14995 0.10314 1.454 0.149
## dlogAM.l1 -0.44881 0.08880 -5.054 1.68e-06 ***
## dlogSEC.l1 0.02674 0.03020 0.885 0.378
## dlogEXT.l1 0.03386 0.05000 0.677 0.500
## const 0.00334 0.01416 0.236 0.814
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.1536 on 113 degrees of freedom
## Multiple R-Squared: 0.1954, Adjusted R-squared: 0.1669
## F-statistic: 6.859 on 4 and 113 DF, p-value: 5.58e-05
##
##
## Estimation results for equation dlogSEC:
## ========================================
## dlogSEC = dlogVM.l1 + dlogAM.l1 + dlogSEC.l1 + dlogEXT.l1 + const
##
## Estimate Std. Error t value Pr(>|t|)
## dlogVM.l1 -0.06904 0.29759 -0.232 0.817
## dlogAM.l1 -0.11140 0.25621 -0.435 0.665
## dlogSEC.l1 -0.45090 0.08714 -5.175 1e-06 ***
## dlogEXT.l1 0.06489 0.14427 0.450 0.654
## const 0.00807 0.04085 0.198 0.844
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.4433 on 113 degrees of freedom
## Multiple R-Squared: 0.2041, Adjusted R-squared: 0.1759
## F-statistic: 7.243 on 4 and 113 DF, p-value: 3.146e-05
##
##
## Estimation results for equation dlogEXT:
## ========================================
## dlogEXT = dlogVM.l1 + dlogAM.l1 + dlogSEC.l1 + dlogEXT.l1 + const
##
## Estimate Std. Error t value Pr(>|t|)
## dlogVM.l1 0.063537 0.195724 0.325 0.7461
## dlogAM.l1 -0.179835 0.168508 -1.067 0.2881
## dlogSEC.l1 0.024311 0.057309 0.424 0.6722
## dlogEXT.l1 -0.168274 0.094888 -1.773 0.0789 .
## const 0.007167 0.026865 0.267 0.7901
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.2916 on 113 degrees of freedom
## Multiple R-Squared: 0.03341, Adjusted R-squared: -0.0008093
## F-statistic: 0.9763 on 4 and 113 DF, p-value: 0.4234
##
##
##
## Covariance matrix of residuals:
## dlogVM dlogAM dlogSEC dlogEXT
## dlogVM 0.019050 0.007519 -0.002253 0.005041
## dlogAM 0.007519 0.023606 0.009789 -0.004818
## dlogSEC -0.002253 0.009789 0.196524 0.017977
## dlogEXT 0.005041 -0.004818 0.017977 0.085011
##
## Correlation matrix of residuals:
## dlogVM dlogAM dlogSEC dlogEXT
## dlogVM 1.00000 0.3545 -0.03682 0.1253
## dlogAM 0.35454 1.0000 0.14372 -0.1075
## dlogSEC -0.03682 0.1437 1.00000 0.1391
## dlogEXT 0.12526 -0.1075 0.13908 1.0000### Se considera 1 rezago#Pruebas del modelo
#Autocorrelación
seriala<-serial.test(Var1, lags.pt=1, type="PT.asymptotic")
seriala$serial ##
## Portmanteau Test (asymptotic)
##
## data: Residuals of VAR object Var1
## Chi-squared = 7.665, df = 0, p-value < 2.2e-16## Posiblemente no hay correlación en los residuos
####Prueba de normalidad
normalidad<-normality.test(Var1)
normalidad$jb.mul## $JB
##
## JB-Test (multivariate)
##
## data: Residuals of VAR object Var1
## Chi-squared = 771.62, df = 8, p-value < 2.2e-16
##
##
## $Skewness
##
## Skewness only (multivariate)
##
## data: Residuals of VAR object Var1
## Chi-squared = 72.676, df = 4, p-value = 6.217e-15
##
##
## $Kurtosis
##
## Kurtosis only (multivariate)
##
## data: Residuals of VAR object Var1
## Chi-squared = 698.94, df = 4, p-value < 2.2e-16## Nohay normalidad en los residuales
#Para la prueba de heteroscedasticidad
arch1<-arch.test(Var1, lags.multi = 8)
arch1$arch.mul##
## ARCH (multivariate)
##
## data: Residuals of VAR object Var1
## Chi-squared = 835.06, df = 800, p-value = 0.1893## Hay heterocedasticidad en los residuales
###Pronóstico modelo VAR
var1.prd <- predict(Var1, n.ahead = 12, ci = 0.95)
#Austar margenes para visualizar
par(mar = c(4, 4, 2, 2))
fanchart(var1.prd)
#Vidrio_manufactura
fanchart(var1.prd, names = "dlogVM",
main = "Fanchart: Predicción de dlogVM",
xlab = "Meses hacia adelante",
ylab = "log Vidrio_manufactura")
#Aluminio_manufacturas
fanchart(var1.prd, names = "dlogAM",
main = "Fanchart: Predicción de dlogAM ",
xlab = "Meses hacia adelante",
ylab = "log Aluminio_manufacturas")
#SECUESTRO
fanchart(var1.prd, names = "dlogSEC",
main = "Fanchart: Predicción de dlogSEC",
xlab = "Meses hacia adelante",
ylab = "log SECUESTRO")
#EXTORSION
fanchart(var1.prd, names = "dlogEXT",
main = "Fanchart: Predicción de dlogEXT",
xlab = "Meses hacia adelante",
ylab = "log EXTORSION")
#Extraer predicciones de dlogVM
dlogVM_pred <- var1.prd$fcst$dlogVM[, "fcst"]
#Último valor real observado de dlogVM
ultimo_dlogVM <- as.numeric(tail(logVM_ts , 1))
#Reconstruir dlogVM proyectado
dlogVM_pred <- cumsum(c(ultimo_dlogVM, dlogVM_pred ))[-1]
#Transformar a escala original
VM_pre <- exp(dlogVM_pred)
#Meses
meses_es <- c("Enero", "Febrero", "Marzo", "Abril", "Mayo", "Junio",
"Julio", "Agosto", "Septiembre", "Octubre", "Noviembre", "Diciembre")
#Crear tabla con resultados
tabla_predicciones <- data.frame(
Mes = meses_es,
Vidrio_manufactura_Proyectada = round(VM_pre, 2)
)
#Ver tabla
print(tabla_predicciones)## Mes Vidrio_manufactura_Proyectada
## 1 Enero 26.06
## 2 Febrero 26.60
## 3 Marzo 26.38
## 4 Abril 26.57
## 5 Mayo 26.55
## 6 Junio 26.63
## 7 Julio 26.67
## 8 Agosto 26.73
## 9 Septiembre 26.78
## 10 Octubre 26.83
## 11 Noviembre 26.88
## 12 Diciembre 26.93### METODOLOGÍA ENGLE-GRANGER para mis cuatro variables
### aplicamos logaritmos a las variables
logVM<-log(Vidrio_manufactura)
logAM<-log(Aluminio_manufacturas)
logSEC<-log(SECUESTRO)
logEXT<-log(EXTORSION)
#Se crean las series de tiempo
x<- ts(logVM, st = c (2015,1),freq=12)
y<- ts(logAM, st = c (2015,1),freq=12)
z<- ts(logSEC, st = c (2015,1),freq=12)
w<- ts(logEXT, st = c (2015,1),freq=12)
## Se organiza en una tabla
datos2<-cbind(x,y,z,w)
datos2## x y z w
## Jan 2015 2.811127 3.376122 2.9957323 6.480045
## Feb 2015 3.104795 3.668274 2.7080502 6.476972
## Mar 2015 3.044327 3.540835 2.8903718 6.352629
## Apr 2015 3.034606 3.553503 2.1972246 6.388561
## May 2015 2.969722 3.536636 2.6390573 6.411818
## Jun 2015 3.036847 3.445289 2.8332133 6.308098
## Jul 2015 3.079626 3.569274 2.8903718 6.289716
## Aug 2015 3.008512 3.415312 2.7725887 6.115892
## Sep 2015 2.960982 3.554648 3.0910425 5.834811
## Oct 2015 2.897567 3.583570 2.1972246 5.429346
## Nov 2015 2.852864 3.478009 2.3025851 5.303305
## Dec 2015 2.813021 3.456299 1.7917595 4.477337
## Jan 2016 2.730539 3.194849 2.7725887 6.431331
## Feb 2016 2.802693 3.386146 2.8903718 6.530878
## Mar 2016 2.831607 3.505445 3.0445224 6.327937
## Apr 2016 2.805693 3.377954 3.3672958 6.295266
## May 2016 2.745752 3.380809 2.6390573 5.880533
## Jun 2016 2.881013 3.459534 2.0794415 5.730100
## Jul 2016 2.819696 3.470016 2.6390573 5.752573
## Aug 2016 2.906600 3.400779 2.8903718 5.796058
## Sep 2016 2.933465 3.494103 2.0794415 5.831882
## Oct 2016 2.828699 3.313342 2.6390573 5.493061
## Nov 2016 2.839225 3.605848 2.5649494 5.433722
## Dec 2016 2.863978 3.477869 2.3978953 5.916202
## Jan 2017 2.738723 3.416516 2.5649494 6.200509
## Feb 2017 2.802609 3.341482 2.5649494 6.165418
## Mar 2017 2.920760 3.444177 3.1780538 6.410175
## Apr 2017 2.887885 3.271490 2.7725887 6.202536
## May 2017 2.887477 3.277493 3.1354942 6.371612
## Jun 2017 2.892213 3.589201 2.3025851 6.182085
## Jul 2017 2.882916 3.321503 2.8903718 6.315358
## Aug 2017 2.796845 3.707756 2.6390573 6.311735
## Sep 2017 2.698948 3.533310 2.5649494 6.173786
## Oct 2017 2.804030 3.667954 2.7080502 5.707110
## Nov 2017 2.804194 3.878257 2.1972246 5.652489
## Dec 2017 2.657851 3.644338 2.5649494 5.402677
## Jan 2018 2.694983 3.550949 2.3978953 6.450470
## Feb 2018 2.632335 3.510768 1.9459101 6.429719
## Mar 2018 2.657631 3.580615 2.4849066 6.391917
## Apr 2018 2.768431 3.623886 3.0445224 6.419995
## May 2018 2.828131 3.555705 2.1972246 6.405228
## Jun 2018 2.780191 3.562223 2.4849066 6.329721
## Jul 2018 2.864944 3.684481 2.4849066 6.576470
## Aug 2018 3.062431 3.607660 2.7080502 6.480045
## Sep 2018 2.794779 3.505854 2.5649494 6.456770
## Oct 2018 2.991363 3.816912 2.3025851 6.350886
## Nov 2018 2.906351 3.694817 2.3025851 6.246107
## Dec 2018 2.863717 3.616011 2.8332133 5.774552
## Jan 2019 2.947657 3.477216 2.6390573 6.679599
## Feb 2019 2.731459 3.422895 1.6094379 6.787845
## Mar 2019 2.847246 3.477889 2.5649494 6.643790
## Apr 2019 3.039786 3.476538 2.3025851 6.453625
## May 2019 3.069801 3.445191 2.1972246 6.535241
## Jun 2019 3.003706 3.601516 2.3025851 6.447306
## Jul 2019 3.058665 3.532046 2.3978953 6.466145
## Aug 2019 2.955961 3.540715 2.4849066 6.486161
## Sep 2019 2.916032 3.429827 1.9459101 6.625392
## Oct 2019 3.036250 3.512130 2.4849066 6.568078
## Nov 2019 2.889966 3.596316 2.3978953 6.469250
## Dec 2019 2.748057 3.614707 2.8903718 6.302619
## Jan 2020 2.895848 3.444708 2.1972246 6.411818
## Feb 2020 2.857688 3.406366 2.4849066 6.423247
## Mar 2020 2.527929 3.379287 2.7080502 5.991465
## Apr 2020 2.840734 2.941659 0.6931472 6.102559
## May 2020 2.127612 3.160527 1.7917595 6.196444
## Jun 2020 2.041221 2.946538 2.7080502 6.259581
## Jul 2020 2.512674 3.388944 2.3978953 6.415097
## Aug 2020 2.543125 3.191355 3.0445224 6.445720
## Sep 2020 2.388557 3.291753 2.4849066 6.553933
## Oct 2020 2.829028 3.630756 2.3025851 6.802395
## Nov 2020 2.897321 3.581910 2.8332133 7.064759
## Dec 2020 2.924057 3.747057 2.1972246 6.993933
## Jan 2021 2.679806 3.503153 2.0794415 6.429719
## Feb 2021 2.762120 3.719353 2.0794415 6.388561
## Mar 2021 2.965582 3.714608 1.7917595 6.320768
## Apr 2021 2.947059 3.600301 2.3978953 6.300786
## May 2021 2.956881 3.384916 2.3025851 6.368187
## Jun 2021 3.122483 3.726347 2.9444390 6.466145
## Jul 2021 3.150100 3.703643 2.3025851 6.552508
## Aug 2021 3.188140 3.954048 2.9444390 6.626718
## Sep 2021 3.266177 3.769816 2.0794415 6.735780
## Oct 2021 3.437052 4.047247 1.7917595 6.736967
## Nov 2021 3.535478 3.894766 1.9459101 6.723832
## Dec 2021 3.654292 4.046820 2.1972246 6.717805
## Jan 2022 3.636922 4.125692 2.3978953 6.376727
## Feb 2022 3.478867 4.083901 2.3025851 6.354370
## Mar 2022 3.474364 3.929422 1.9459101 6.487684
## Apr 2022 3.451833 3.793360 1.7917595 6.495266
## May 2022 3.431541 3.907165 2.4849066 6.563856
## Jun 2022 3.654402 4.298532 2.6390573 6.579251
## Jul 2022 3.593020 4.279853 2.3978953 6.536692
## Aug 2022 3.686396 4.301708 2.7725887 6.786717
## Sep 2022 3.666231 4.088317 2.5649494 6.886532
## Oct 2022 3.667945 4.188050 3.0910425 7.000334
## Nov 2022 3.519770 4.140202 3.0910425 6.951772
## Dec 2022 3.485445 3.855891 2.9444390 7.075809
## Jan 2023 3.355175 3.744458 3.2958369 6.563856
## Feb 2023 3.339672 3.752219 2.7725887 6.565265
## Mar 2023 3.261369 4.000926 3.3322045 6.727432
## Apr 2023 3.190696 3.767283 2.8332133 6.669498
## May 2023 3.247684 3.783881 3.0910425 6.747587
## Jun 2023 3.120841 3.781255 2.8332133 6.646391
## Jul 2023 3.008149 3.721235 3.1354942 6.701960
## Aug 2023 3.114891 3.693488 3.2188758 6.875232
## Sep 2023 2.956411 3.558887 2.9444390 6.971669
## Oct 2023 3.072139 3.730474 3.4657359 7.018402
## Nov 2023 3.011003 3.681716 2.8332133 7.135687
## Dec 2023 3.001119 3.825041 2.9444390 7.084226
## Jan 2024 3.147920 3.755495 3.2188758 6.742881
## Feb 2024 3.189283 3.602555 3.0910425 6.865891
## Mar 2024 2.897548 3.586814 3.0910425 6.903747
## Apr 2024 3.229888 3.744570 3.0910425 6.926577
## May 2024 3.106644 3.903858 2.9444390 7.068172
## Jun 2024 3.061325 3.472076 2.8332133 6.966024
## Jul 2024 3.139821 3.728326 3.4339872 7.104965
## Aug 2024 3.260230 3.902465 3.2188758 7.164720
## Sep 2024 3.187717 3.796556 2.8903718 7.158514
## Oct 2024 3.397241 3.896235 3.6109179 7.197435
## Nov 2024 3.283876 3.975313 3.6888795 7.210818
## Dec 2024 3.284544 3.907596 3.2958369 7.138867class(datos2)## [1] "mts" "ts" "matrix" "array"### Se grafican
ts.plot(datos2, col=c("blue","red"), main = "Tendencia")
#Generar Modelo
modelo1<-lm(y~x+z+w)
summary(modelo1)##
## Call:
## lm(formula = y ~ x + z + w)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.52735 -0.11182 -0.01033 0.10767 0.41482
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.36094 0.19737 6.895 2.99e-10 ***
## x 0.67942 0.04889 13.898 < 2e-16 ***
## z 0.02083 0.03063 0.680 0.498
## w 0.02680 0.03468 0.773 0.441
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1439 on 116 degrees of freedom
## Multiple R-squared: 0.6889, Adjusted R-squared: 0.6808
## F-statistic: 85.6 on 3 and 116 DF, p-value: < 2.2e-16residualesmodelo1<-modelo1$residuals
summary(residualesmodelo1)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.52735 -0.11182 -0.01033 0.00000 0.10767 0.41482library(car)
residualPlots(modelo1,fitted=TRUE)
## Test stat Pr(>|Test stat|)
## x 2.5580 0.011827 *
## z -3.0326 0.002997 **
## w 1.5177 0.131842
## Tukey test 2.6619 0.007771 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Evaluando la estacionaridad de los residuales
library(tseries)
adf.test(residualesmodelo1)##
## Augmented Dickey-Fuller Test
##
## data: residualesmodelo1
## Dickey-Fuller = -3.5681, Lag order = 4, p-value = 0.03913
## alternative hypothesis: stationary#Rechazo H0 por lo tanto los residuales son estacionarios
#las variables estan cointegrada
#Con tendencia
library(urca)
adf.test1<-ur.df(residualesmodelo1,type="trend",selectlags = "AIC")
summary(adf.test1)##
## ###############################################
## # 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.51990 -0.08708 -0.00109 0.08682 0.36104
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0260016 0.0257068 -1.011 0.3139
## z.lag.1 -0.5823399 0.1117287 -5.212 8.41e-07 ***
## tt 0.0004557 0.0003725 1.223 0.2237
## z.diff.lag -0.2215558 0.0912533 -2.428 0.0167 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1341 on 114 degrees of freedom
## Multiple R-squared: 0.4049, Adjusted R-squared: 0.3893
## F-statistic: 25.86 on 3 and 114 DF, p-value: 7.815e-13
##
##
## Value of test-statistic is: -5.2121 9.0624 13.583
##
## 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.47adf.test1@teststat## tau3 phi2 phi3
## statistic -5.212088 9.062435 13.58304adf.test1@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### -5.2121 es menor -3.43, no esta zona de rechazo por tanto rechazo H0, la serie es estacionaria
#con constante
adf.test2<-ur.df(residualesmodelo1,type="drift",selectlags = "AIC")
summary(adf.test2)##
## ###############################################
## # 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.51817 -0.08069 -0.00971 0.09416 0.34401
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.001584 0.012371 0.128 0.89836
## z.lag.1 -0.550630 0.108915 -5.056 1.64e-06 ***
## z.diff.lag -0.237524 0.090510 -2.624 0.00986 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1344 on 115 degrees of freedom
## Multiple R-squared: 0.3971, Adjusted R-squared: 0.3866
## F-statistic: 37.87 on 2 and 115 DF, p-value: 2.31e-13
##
##
## Value of test-statistic is: -5.0556 12.7901
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1 6.52 4.63 3.81#-5.0556 menor a -2.88, cae en zona de no rechazo por tanto rechazo H0
#entonces la serie es estacionaria
#con constante y tendencia
adf.test3<-ur.df(residualesmodelo1,type="none",selectlags = "AIC")
summary(adf.test3)##
## ###############################################
## # 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.51661 -0.07911 -0.00813 0.09576 0.34559
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 -0.55068 0.10845 -5.078 1.47e-06 ***
## z.diff.lag -0.23739 0.09012 -2.634 0.00959 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1338 on 116 degrees of freedom
## Multiple R-squared: 0.397, Adjusted R-squared: 0.3866
## F-statistic: 38.19 on 2 and 116 DF, p-value: 1.807e-13
##
##
## Value of test-statistic is: -5.0776
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.58 -1.95 -1.62#-5.0776 menor a -1.95, cae en zona de no rechazo por tanto rechazo H0
#entonces la serie es estacionaria
#De acuerdo con esta revisión, los residuos del modelo muestran que hay cointegración, situación que beneficia la elaboración del modelo VEC.
# Prueba de phillips - ouliaris para cointegración
#Con p2 no importa cual variable esta al lado izquierdo de la ecuación
prueba.PO<-ca.po(datos2,type="Pz")
summary(prueba.PO)##
## ########################################
## # Phillips and Ouliaris Unit Root Test #
## ########################################
##
## Test of type Pz
## detrending of series none
##
## Response x :
##
## Call:
## lm(formula = x ~ zr - 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.55879 -0.07314 -0.00255 0.07456 0.32005
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## zrx 0.70431 0.07344 9.590 2.33e-16 ***
## zry 0.18671 0.07507 2.487 0.0143 *
## zrz 0.03769 0.02955 1.275 0.2047
## zrw 0.01807 0.02802 0.645 0.5203
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.138 on 115 degrees of freedom
## Multiple R-squared: 0.998, Adjusted R-squared: 0.9979
## F-statistic: 1.419e+04 on 4 and 115 DF, p-value: < 2.2e-16
##
##
## Response y :
##
## Call:
## lm(formula = y ~ zr - 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.42311 -0.09428 -0.01077 0.09208 0.39608
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## zrx 0.22999 0.08113 2.835 0.00542 **
## zry 0.57244 0.08293 6.903 2.96e-10 ***
## zrz -0.01055 0.03264 -0.323 0.74705
## zrw 0.13823 0.03096 4.465 1.88e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1525 on 115 degrees of freedom
## Multiple R-squared: 0.9983, Adjusted R-squared: 0.9982
## F-statistic: 1.691e+04 on 4 and 115 DF, p-value: < 2.2e-16
##
##
## Response z :
##
## Call:
## lm(formula = z ~ zr - 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.7805 -0.2361 0.0519 0.2766 0.9859
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## zrx 0.22455 0.21765 1.032 0.30438
## zry -0.19476 0.22247 -0.875 0.38315
## zrz 0.37923 0.08757 4.331 3.19e-05 ***
## zrw 0.25655 0.08305 3.089 0.00252 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.409 on 115 degrees of freedom
## Multiple R-squared: 0.977, Adjusted R-squared: 0.9762
## F-statistic: 1221 on 4 and 115 DF, p-value: < 2.2e-16
##
##
## Response w :
##
## Call:
## lm(formula = w ~ zr - 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.9622 -0.1079 0.0195 0.1024 1.7046
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## zrx -0.12305 0.15090 -0.815 0.417
## zry 0.36083 0.15424 2.339 0.021 *
## zrz 0.02312 0.06071 0.381 0.704
## zrw 0.84520 0.05758 14.679 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2836 on 115 degrees of freedom
## Multiple R-squared: 0.9981, Adjusted R-squared: 0.9981
## F-statistic: 1.541e+04 on 4 and 115 DF, p-value: < 2.2e-16
##
##
##
## Value of test-statistic is: 280.5945
##
## Critical values of Pz are:
## 10pct 5pct 1pct
## critical values 99.2664 109.7426 131.5716#Si el test estadistico es menor al valor critico 5%
# las variables test-statistic=280.5945 mayor al valor critico=109.7426
#Por lo tanto las variables estan cointegradas
prueba.PU<-ca.po(datos2,type="Pu")
summary(prueba.PU)##
## ########################################
## # 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
## -0.50877 -0.11300 0.02451 0.11711 0.36458
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z[, -1]y 0.80667 0.05794 13.922 <2e-16 ***
## z[, -1]z -0.01646 0.03703 -0.445 0.657
## z[, -1]w 0.01878 0.03525 0.533 0.595
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1738 on 117 degrees of freedom
## Multiple R-squared: 0.9968, Adjusted R-squared: 0.9967
## F-statistic: 1.2e+04 on 3 and 117 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: 64.8923
##
## Critical values of Pu are:
## 10pct 5pct 1pct
## critical values 33.5359 40.122 55.7341#Si el test estadistico es menor al valor critico 5%
# las variables test-statistic= 64.8923 mayor al valor critico=40.122
#Por lo tanto las variables estan cointegradas
###################### MODELO DE CORRECCION DE ERRORES ###################
#Segundo paso: Generar el modelo de correción de errores
#Generar la variable tendencia para agregarla al modelo y revisar si
#genera residuales estacionariostendencia <- seq_along(y)
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 120modelo2<-lm(y~tendencia+x+z+w)
summary(modelo2)##
## Call:
## lm(formula = y ~ tendencia + x + z + w)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.56677 -0.07805 -0.01011 0.08515 0.38214
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.9700166 0.2616635 7.529 1.25e-11 ***
## tendencia 0.0018522 0.0005499 3.368 0.00103 **
## x 0.6300414 0.0490842 12.836 < 2e-16 ***
## z 0.0238839 0.0293681 0.813 0.41775
## w -0.0633282 0.0426682 -1.484 0.14049
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1379 on 115 degrees of freedom
## Multiple R-squared: 0.7168, Adjusted R-squared: 0.7069
## F-statistic: 72.76 on 4 and 115 DF, p-value: < 2.2e-16modelo2##
## Call:
## lm(formula = y ~ tendencia + x + z + w)
##
## Coefficients:
## (Intercept) tendencia x z w
## 1.970017 0.001852 0.630041 0.023884 -0.063328residualesmodelo2 <- modelo2$residuals
summary(residualesmodelo2)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.56677 -0.07805 -0.01011 0.00000 0.08515 0.38214residualPlots(modelo2)
## Test stat Pr(>|Test stat|)
## tendencia -1.7825 0.07733 .
## x 1.8220 0.07108 .
## z -4.0725 8.62e-05 ***
## w 0.2770 0.78229
## Tukey test 1.9527 0.05085 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1modelo2##
## Call:
## lm(formula = y ~ tendencia + x + z + w)
##
## Coefficients:
## (Intercept) tendencia x z w
## 1.970017 0.001852 0.630041 0.023884 -0.063328# Probar que los residuales son estacionarios
# Prueba Dickey-fuller para revisar estacionariedad
adf.test2 <- ur.df(residualesmodelo2, type="none",selectlags = "AIC")
summary(adf.test2)##
## ###############################################
## # 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.53913 -0.08234 0.00159 0.09250 0.31872
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 -0.56842 0.10990 -5.172 9.78e-07 ***
## z.diff.lag -0.23233 0.09025 -2.574 0.0113 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1287 on 116 degrees of freedom
## Multiple R-squared: 0.4052, Adjusted R-squared: 0.395
## F-statistic: 39.52 on 2 and 116 DF, p-value: 8.164e-14
##
##
## Value of test-statistic is: -5.1723
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.58 -1.95 -1.62##-5.1723 es menor que -1.95 Acepto H0:
# Los residuales son estacionarios al 1% 5% y 10%
##Cuantas diferencias se necesitan
library(forecast)
ndiffs(x)#1## [1] 1ndiffs(y)#1## [1] 1ndiffs(z)#1## [1] 1ndiffs(w)#1## [1] 1#Para calcular el modelo de correción de errores se generan las diferencias
# y el error con un rezago
dy<-diff(y)
dx<-diff(x)
dz<-diff(z)
dw<-diff(w)
#si una variable no es significativa se elimina
modelo3<-lm(dy~dx+dz+dw)
summary(modelo3)##
## Call:
## lm(formula = dy ~ dx + dz + dw)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.42093 -0.11881 -0.00367 0.09545 0.42990
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.00327 0.01481 0.221 0.82559
## dx 0.36564 0.10286 3.555 0.00055 ***
## dz 0.06253 0.03112 2.010 0.04682 *
## dw -0.07516 0.05206 -1.444 0.15150
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1614 on 115 degrees of freedom
## Multiple R-squared: 0.1183, Adjusted R-squared: 0.09529
## F-statistic: 5.143 on 3 and 115 DF, p-value: 0.002255###eliminamos w=extorsion que no es significativa para el modelo
modelo31 <- lm(dy ~ dx + dz)
summary(modelo31)##
## Call:
## lm(formula = dy ~ dx + dz)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.43497 -0.11808 -0.00285 0.09797 0.42719
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.00295 0.01487 0.198 0.84310
## dx 0.34527 0.10236 3.373 0.00101 **
## dz 0.05639 0.03097 1.821 0.07121 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1622 on 116 degrees of freedom
## Multiple R-squared: 0.1023, Adjusted R-squared: 0.08683
## F-statistic: 6.61 on 2 and 116 DF, p-value: 0.001912anova(modelo31, modelo3)## Analysis of Variance Table
##
## Model 1: dy ~ dx + dz
## Model 2: dy ~ dx + dz + dw
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 116 3.0510
## 2 115 2.9967 1 0.054323 2.0847 0.1515#acepto H0, La variable dw no contribuye significativamente a explicar la variable dependiente dy.
#Por tanto, eliminamos dw del modelo y trabajar con el modelo reducido dy ~ dx + dz.
#Se generan los residuales
res3<-modelo31$residuals
# Se calcula un rezago en los residuales
res3_1 <- lag(res3) #Un rezago
# Se genera el modelo de correcion de errores
modelo4<-lm(dy~dx+res3_1)
summary(modelo4)##
## Call:
## lm(formula = dy ~ dx + res3_1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.51279 -0.08588 0.00272 0.08407 0.35093
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.001703 0.013260 0.128 0.89806
## dx 0.345414 0.091910 3.758 0.00027 ***
## res3_1 -0.485528 0.082845 -5.861 4.49e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.144 on 115 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.2803, Adjusted R-squared: 0.2678
## F-statistic: 22.4 on 2 and 115 DF, p-value: 6.089e-09MCE1 <- modelo4
summary(MCE1)##
## Call:
## lm(formula = dy ~ dx + res3_1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.51279 -0.08588 0.00272 0.08407 0.35093
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.001703 0.013260 0.128 0.89806
## dx 0.345414 0.091910 3.758 0.00027 ***
## res3_1 -0.485528 0.082845 -5.861 4.49e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.144 on 115 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.2803, Adjusted R-squared: 0.2678
## F-statistic: 22.4 on 2 and 115 DF, p-value: 6.089e-09## Veamos la estacionariedad de las series:
library(tseries)
#Prueba Phillips-Perron
pp.test(y)## Warning in pp.test(y): p-value smaller than printed p-value##
## Phillips-Perron Unit Root Test
##
## data: y
## Dickey-Fuller Z(alpha) = -36.751, Truncation lag parameter = 4, p-value
## = 0.01
## alternative hypothesis: stationary# p-value = 0.01 La serie es estacionaria
pp.test(x)##
## Phillips-Perron Unit Root Test
##
## data: x
## Dickey-Fuller Z(alpha) = -17.169, Truncation lag parameter = 4, p-value
## = 0.1228
## alternative hypothesis: stationary# p-value = 0.1228La serie es estacionaria
#Veamos la cointegración
#Test Phillips-Ouliaris Cointegration
po.test(cbind(x,y,z,w))## Warning in po.test(cbind(x, y, z, w)): p-value smaller than printed p-value##
## Phillips-Ouliaris Cointegration Test
##
## data: cbind(x, y, z, w)
## Phillips-Ouliaris demeaned = -62.526, Truncation lag parameter = 1,
## p-value = 0.01#H0: No Cointegradas
#p-value = 0.01 acepto H0 por tanto las variables estan cointegradas
#Primer paso estimar el modelo:
m1 <- lm(y~x+z+w)
summary(m1)##
## Call:
## lm(formula = y ~ x + z + w)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.52735 -0.11182 -0.01033 0.10767 0.41482
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.36094 0.19737 6.895 2.99e-10 ***
## x 0.67942 0.04889 13.898 < 2e-16 ***
## z 0.02083 0.03063 0.680 0.498
## w 0.02680 0.03468 0.773 0.441
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1439 on 116 degrees of freedom
## Multiple R-squared: 0.6889, Adjusted R-squared: 0.6808
## F-statistic: 85.6 on 3 and 116 DF, p-value: < 2.2e-16# Segundo paso: determinar si los residuas del modelo son I(0)
resm1<-resid(m1)
library(tseries)
pp.test(resm1)## Warning in pp.test(resm1): p-value smaller than printed p-value##
## Phillips-Perron Unit Root Test
##
## data: resm1
## Dickey-Fuller Z(alpha) = -101.44, Truncation lag parameter = 4, p-value
## = 0.01
## alternative hypothesis: stationaryBox.test(resm1,type = "Ljung-Box", lag=5) #Ho: independencia##
## Box-Ljung test
##
## data: resm1
## X-squared = 28.273, df = 5, p-value = 3.219e-05# Ahora estimamos el modelo
n <- length(y)
MCE2<-lm(diff(y)~diff(x)+lag(resm1[1:(n-1)],1)) #se incluyen 5 rezagos
summary(MCE2)##
## Call:
## lm(formula = diff(y) ~ diff(x) + lag(resm1[1:(n - 1)], 1))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.52742 -0.10899 -0.00369 0.10628 0.42490
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.001628 0.015105 0.108 0.91434
## diff(x) 0.295327 0.104260 2.833 0.00545 **
## lag(resm1[1:(n - 1)], 1) 0.035084 0.106211 0.330 0.74175
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1641 on 115 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.06629, Adjusted R-squared: 0.05005
## F-statistic: 4.082 on 2 and 115 DF, p-value: 0.01937Box.test(resid(MCE2)) #hay independencia##
## Box-Pierce test
##
## data: resid(MCE2)
## X-squared = 23.847, df = 1, p-value = 1.043e-06# Se ha corregido la estructura
## Tenemos el efecto a largo plazo
summary(m1)##
## Call:
## lm(formula = y ~ x + z + w)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.52735 -0.11182 -0.01033 0.10767 0.41482
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.36094 0.19737 6.895 2.99e-10 ***
## x 0.67942 0.04889 13.898 < 2e-16 ***
## z 0.02083 0.03063 0.680 0.498
## w 0.02680 0.03468 0.773 0.441
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1439 on 116 degrees of freedom
## Multiple R-squared: 0.6889, Adjusted R-squared: 0.6808
## F-statistic: 85.6 on 3 and 116 DF, p-value: < 2.2e-16summary(MCE1)##
## Call:
## lm(formula = dy ~ dx + res3_1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.51279 -0.08588 0.00272 0.08407 0.35093
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.001703 0.013260 0.128 0.89806
## dx 0.345414 0.091910 3.758 0.00027 ***
## res3_1 -0.485528 0.082845 -5.861 4.49e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.144 on 115 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.2803, Adjusted R-squared: 0.2678
## F-statistic: 22.4 on 2 and 115 DF, p-value: 6.089e-09summary(MCE2)##
## Call:
## lm(formula = diff(y) ~ diff(x) + lag(resm1[1:(n - 1)], 1))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.52742 -0.10899 -0.00369 0.10628 0.42490
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.001628 0.015105 0.108 0.91434
## diff(x) 0.295327 0.104260 2.833 0.00545 **
## lag(resm1[1:(n - 1)], 1) 0.035084 0.106211 0.330 0.74175
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1641 on 115 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.06629, Adjusted R-squared: 0.05005
## F-statistic: 4.082 on 2 and 115 DF, p-value: 0.01937Pronostico
# Cargar librerías necesarias
library(forecast)
library(dplyr)
# Datos originales (suponemos que 'datos' es el data frame completo)
# Crear series transformadas como en el taller
datos$logVM <- log(datos$Vidrio_manufactura / 1000)
datos$logAM <- log(datos$Aluminio_manufacturas / 1000)
datos$logSEC <- log(datos$SECUESTRO)
datos$logEXT <- log(datos$EXTORSION)
# Crear series de tiempo (frecuencia mensual, como en el taller)
logVM_ts <- ts(datos$logVM, start = c(2015, 1), frequency = 12)
logAM_ts <- ts(datos$logAM, start = c(2015, 1), frequency = 12)
logSEC_ts <- ts(datos$logSEC, start = c(2015, 1), frequency = 12)
logEXT_ts <- ts(datos$logEXT, start = c(2015, 1), frequency = 12)
# Calcular primeras diferencias
dx <- diff(logVM_ts)
dy <- diff(logAM_ts)
dz <- diff(logSEC_ts)
dw <- diff(logEXT_ts)
# Crear datos2 como en el taller
datos2 <- cbind(x = logVM_ts, y = logAM_ts, z = logSEC_ts, w = logEXT_ts)
# Paso 1: Pronosticar logAM con MCE1
# Coeficientes de MCE1: dy = 0.001703 + 0.345414 * dx - 0.485528 * res3_1
n_ahead <- 12
dy_pred <- numeric(n_ahead)
ultimo_logAM <- as.numeric(tail(datos2[, "y"], 1))
ultimo_dx <- as.numeric(tail(dx, 1))
ultimo_dz <- as.numeric(tail(dz, 1))
# Ajustar ARIMA para dx y dz
arima_dx <- Arima(dx, order = c(1, 0, 0))
dx_pred <- forecast(arima_dx, h = n_ahead)$mean
arima_dz <- Arima(dz, order = c(1, 0, 0))
dz_pred <- forecast(arima_dz, h = n_ahead)$mean
# Reestimar modelo reducido para res3_1 (dy ~ dx + dz)
modelo_reducido <- lm(dy ~ dx + dz)
beta0 <- coef(modelo_reducido)[1]
beta1 <- coef(modelo_reducido)[2]
beta2 <- coef(modelo_reducido)[3]
# Calcular res3_1 inicial
ultimo_dy <- as.numeric(tail(dy, 1))
ultimo_res3_1 <- ultimo_dy - (beta0 + beta1 * ultimo_dx + beta2 * ultimo_dz)
# Pronosticar dy y res3_1 iterativamente
res3_1_pred <- numeric(n_ahead)
for (i in 1:n_ahead) {
dy_pred[i] <- 0.001703 + 0.345414 * dx_pred[i] - 0.485528 * ultimo_res3_1
res3_1_pred[i] <- dy_pred[i] - (beta0 + beta1 * dx_pred[i] + beta2 * dz_pred[i])
ultimo_res3_1 <- res3_1_pred[i]
}
# Reconstruir logAM
logAM_pred <- ultimo_logAM + cumsum(c(0, dy_pred))[1:n_ahead]
AM_pre <- exp(logAM_pred) * 1000 # Escala original
# Paso 2: Pronosticar logVM, logSEC, y logEXT con ARIMA
arima_dw <- Arima(dw, order = c(1, 0, 0))
dw_pred <- forecast(arima_dw, h = n_ahead)$mean
# Reconstruir series en logaritmos
ultimo_logVM <- as.numeric(tail(datos2[, "x"], 1))
ultimo_logSEC <- as.numeric(tail(datos2[, "z"], 1))
ultimo_logEXT <- as.numeric(tail(datos2[, "w"], 1))
logVM_pred <- ultimo_logVM + cumsum(c(0, dx_pred))[1:n_ahead]
logSEC_pred <- ultimo_logSEC + cumsum(c(0, dz_pred))[1:n_ahead]
logEXT_pred <- ultimo_logEXT + cumsum(c(0, dw_pred))[1:n_ahead]
# Volver a la escala original
VM_pre <- exp(logVM_pred) * 1000
SEC_pre <- exp(logSEC_pred)
EXT_pre <- exp(logEXT_pred)
# Paso 3: Crear tabla de predicciones
meses_es <- c("Enero", "Febrero", "Marzo", "Abril", "Mayo", "Junio",
"Julio", "Agosto", "Septiembre", "Octubre", "Noviembre", "Diciembre")
tabla_predicciones <- data.frame(
Mes = meses_es,
Vidrio_manufactura_Proyectada = round(VM_pre, 2),
Aluminio_manufacturas_Proyectada = round(AM_pre, 2),
SECUESTRO_Proyectado = round(SEC_pre, 2),
EXTORSION_Proyectado = round(EXT_pre, 2)
)
# Mostrar tabla
print(tabla_predicciones)## Mes Vidrio_manufactura_Proyectada Aluminio_manufacturas_Proyectada
## 1 Enero 26.70 49.78
## 2 Febrero 26.81 51.13
## 3 Marzo 26.90 50.97
## 4 Abril 26.99 51.20
## 5 Mayo 27.08 51.39
## 6 Junio 27.18 51.53
## 7 Julio 27.27 51.74
## 8 Agosto 27.36 51.89
## 9 Septiembre 27.46 52.08
## 10 Octubre 27.55 52.25
## 11 Noviembre 27.65 52.43
## 12 Diciembre 27.74 52.60
## SECUESTRO_Proyectado EXTORSION_Proyectado
## 1 27.00 1260.00
## 2 32.37 1282.05
## 3 30.04 1287.01
## 4 31.25 1294.63
## 5 30.90 1301.89
## 6 31.25 1309.25
## 7 31.29 1316.64
## 8 31.47 1324.08
## 9 31.58 1331.56
## 10 31.73 1339.08
## 11 31.86 1346.65
## 12 32.00 1354.25# Validar rango de los pronósticos
cat("Rango de Vidrio_manufactura original:", range(datos$Vidrio_manufactura), "\n")## Rango de Vidrio_manufactura original: 7.700008 39.90077cat("Rango de Aluminio_manufacturas original:", range(datos$Aluminio_manufacturas), "\n")## Rango de Aluminio_manufacturas original: 18.94725 73.8258cat("Rango de SECUESTRO original:", range(datos$SECUESTRO), "\n")## Rango de SECUESTRO original: 2 40cat("Rango de EXTORSION original:", range(datos$EXTORSION), "\n")## Rango de EXTORSION original: 88 1354# Visualizar pronósticos
par(mfrow = c(2, 2))
plot(1:n_ahead, VM_pre, type = "l", col = "blue", lwd = 2,
main = "Pronóstico: Vidrio_manufactura", xlab = "Meses hacia adelante",
ylab = "Vidrio_manufactura", ylim = range(datos$Vidrio_manufactura))
grid()
plot(1:n_ahead, AM_pre, type = "l", col = "blue", lwd = 2,
main = "Pronóstico: Aluminio_manufacturas (MCE1)", xlab = "Meses hacia adelante",
ylab = "Aluminio_manufacturas", ylim = range(datos$Aluminio_manufacturas))
grid()
plot(1:n_ahead, SEC_pre, type = "l", col = "blue", lwd = 2,
main = "Pronóstico: SECUESTRO", xlab = "Meses hacia adelante",
ylab = "SECUESTRO", ylim = range(datos$SECUESTRO))
grid()
plot(1:n_ahead, EXT_pre, type = "l", col = "blue", lwd = 2,
main = "Pronóstico: EXTORSION", xlab = "Meses hacia adelante",
ylab = "EXTORSION", ylim = range(datos$EXTORSION))
grid()
METODOLOGÍA JOHANSEN-JUSELIOUS
logVM<-log(Vidrio_manufactura)
logAM<-log(Aluminio_manufacturas)
logSEC<-log(SECUESTRO)
logEXT<-log(EXTORSION)
datos3 <- cbind(logVM,logAM,logSEC,logEXT)# Calculando el numero optimo de retardos
library(vars)
VARselect(datos3, lag.max = 12)## $selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 1 1 1 1
##
## $criteria
## 1 2 3 4 5
## AIC(n) -1.225514e+01 -1.216935e+01 -1.210715e+01 -1.206658e+01 -1.193489e+01
## HQ(n) -1.205375e+01 -1.180685e+01 -1.158353e+01 -1.138185e+01 -1.108906e+01
## SC(n) -1.175845e+01 -1.127531e+01 -1.081575e+01 -1.037783e+01 -9.848793e+00
## FPE(n) 4.761846e-06 5.195072e-06 5.545884e-06 5.809468e-06 6.690425e-06
## 6 7 8 9 10
## AIC(n) -1.183891e+01 -1.179638e+01 -1.171739e+01 -1.171401e+01 -1.171341e+01
## HQ(n) -1.083196e+01 -1.062832e+01 -1.038821e+01 -1.022372e+01 -1.006201e+01
## SC(n) -9.355455e+00 -8.915576e+00 -8.439226e+00 -8.038496e+00 -7.640546e+00
## FPE(n) 7.469166e-06 7.949518e-06 8.835403e-06 9.179527e-06 9.603006e-06
## 11 12
## AIC(n) -1.161363e+01 -1.170024e+01
## HQ(n) -9.801122e+00 -9.726619e+00
## SC(n) -7.143415e+00 -6.832671e+00
## FPE(n) 1.122065e-05 1.102838e-05# En primer lugar se aplica la prueba de la traza de cointegración
# con tendencia y constante
# Calculando el numero optimo de retardos
library(vars)
VARselect(datos3, lag.max = 12)## $selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 1 1 1 1
##
## $criteria
## 1 2 3 4 5
## AIC(n) -1.225514e+01 -1.216935e+01 -1.210715e+01 -1.206658e+01 -1.193489e+01
## HQ(n) -1.205375e+01 -1.180685e+01 -1.158353e+01 -1.138185e+01 -1.108906e+01
## SC(n) -1.175845e+01 -1.127531e+01 -1.081575e+01 -1.037783e+01 -9.848793e+00
## FPE(n) 4.761846e-06 5.195072e-06 5.545884e-06 5.809468e-06 6.690425e-06
## 6 7 8 9 10
## AIC(n) -1.183891e+01 -1.179638e+01 -1.171739e+01 -1.171401e+01 -1.171341e+01
## HQ(n) -1.083196e+01 -1.062832e+01 -1.038821e+01 -1.022372e+01 -1.006201e+01
## SC(n) -9.355455e+00 -8.915576e+00 -8.439226e+00 -8.038496e+00 -7.640546e+00
## FPE(n) 7.469166e-06 7.949518e-06 8.835403e-06 9.179527e-06 9.603006e-06
## 11 12
## AIC(n) -1.161363e+01 -1.170024e+01
## HQ(n) -9.801122e+00 -9.726619e+00
## SC(n) -7.143415e+00 -6.832671e+00
## FPE(n) 1.122065e-05 1.102838e-05# En primer lugar se aplica la prueba de la traza de cointegración
# con tendencia y constante
library(urca)
summary(ca.jo(datos3,type = "trace",ecdet = "none",spec = c("longrun"),K=2))##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: trace statistic , with linear trend
##
## Eigenvalues (lambda):
## [1] 0.23778663 0.16053132 0.11196139 0.02950211
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 3 | 3.53 6.50 8.18 11.65
## r <= 2 | 17.54 15.66 17.95 23.52
## r <= 1 | 38.19 28.71 31.52 37.22
## r = 0 | 70.23 45.23 48.28 55.43
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## logVM.l2 logAM.l2 logSEC.l2 logEXT.l2
## logVM.l2 1.0000000 1.0000000 1.00000000 1.00000000
## logAM.l2 -1.5337623 -1.1270718 0.08946243 0.77995849
## logSEC.l2 -0.4091275 0.3420416 -0.61904255 0.09143807
## logEXT.l2 0.4208877 -0.2642130 -0.79299100 0.28933033
##
## Weights W:
## (This is the loading matrix)
##
## logVM.l2 logAM.l2 logSEC.l2 logEXT.l2
## logVM.d -0.04576371 -0.169638120 -0.05456330 -0.02416623
## logAM.d 0.19703225 0.061873018 -0.04720009 -0.02984947
## logSEC.d 0.57948408 -0.451084842 0.18452091 -0.01718318
## logEXT.d -0.26145793 0.001660205 0.12765052 -0.06088780##Hay evidencia de al menos 2 vectores de cointegración al 5% (y posiblemente 3 al 10%).
#Prueba de JOHANSEN de la traza con constante
summary(ca.jo(datos3,type = "trace",ecdet = "const",spec = c("longrun"),K=2))##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: trace statistic , without linear trend and constant in cointegration
##
## Eigenvalues (lambda):
## [1] 2.378559e-01 1.605493e-01 1.120584e-01 3.064810e-02 -1.270190e-18
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 3 | 3.67 7.52 9.24 12.97
## r <= 2 | 17.70 17.85 19.96 24.60
## r <= 1 | 38.35 32.00 34.91 41.07
## r = 0 | 70.40 49.65 53.12 60.16
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## logVM.l2 logAM.l2 logSEC.l2 logEXT.l2 constant
## logVM.l2 1.0000000 1.0000000 1.00000000 1.00000000 1.0000000
## logAM.l2 -1.5340122 -1.1273346 0.08851884 0.76475304 1.2379625
## logSEC.l2 -0.4095895 0.3419297 -0.61826718 0.07754885 0.4730805
## logEXT.l2 0.4209455 -0.2636368 -0.79586595 0.28412769 0.4021763
## constant 0.9224470 1.8922082 3.42289149 -7.90692840 -7.6999659
##
## Weights W:
## (This is the loading matrix)
##
## logVM.l2 logAM.l2 logSEC.l2 logEXT.l2 constant
## logVM.d -0.04552778 -0.169768303 -0.05429868 -0.02453659 1.113079e-16
## logAM.d 0.19713086 0.061864066 -0.04687007 -0.03026915 -6.601442e-17
## logSEC.d 0.57999898 -0.451568236 0.18421321 -0.01690699 1.079698e-16
## logEXT.d -0.26049455 0.000765155 0.12816694 -0.06147254 -1.150006e-17##Hay evidencia de 2 vectores de cointegración al 5%.
#Prueba de JOHANSEN de la traza con tendencia
summary(ca.jo(datos3,type = "trace",ecdet = "trend",spec = c("longrun"),K=2))##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: trace statistic , with linear trend in cointegration
##
## Eigenvalues (lambda):
## [1] 2.713727e-01 2.201703e-01 1.466493e-01 6.366977e-02 1.665335e-16
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 3 | 7.76 10.49 12.25 16.26
## r <= 2 | 26.48 22.76 25.32 30.45
## r <= 1 | 55.82 39.06 42.44 48.45
## r = 0 | 93.18 59.14 62.99 70.05
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## logVM.l2 logAM.l2 logSEC.l2 logEXT.l2 trend.l2
## logVM.l2 1.000000 1.000000000 1.000000000 1.000000000 1.00000000
## logAM.l2 -13.018481 -1.408017867 -1.291809322 0.460361002 0.20679965
## logSEC.l2 -65.769764 -0.173870988 0.801990168 -0.122778996 0.48830609
## logEXT.l2 197.088088 -0.133256767 0.078772531 -0.079849817 0.15419215
## trend.l2 -1.609663 0.003369843 -0.003784279 -0.006650682 0.03930157
##
## Weights W:
## (This is the loading matrix)
##
## logVM.l2 logAM.l2 logSEC.l2 logEXT.l2 trend.l2
## logVM.d 0.0002490604 -0.17268951 -0.04094915 -0.08554622 -7.741429e-17
## logAM.d 0.0004675289 0.23715128 0.03782472 -0.08490209 -8.010488e-17
## logSEC.d 0.0013239062 0.66119039 -0.37015577 0.02774634 -3.414758e-16
## logEXT.d -0.0024808610 0.09894332 -0.13771592 -0.05503824 5.661135e-16##Hay evidencia de 3 vectores de cointegración al 5%.
#Con type =eigen
# con tendencia y constante
library(urca)
summary(ca.jo(datos3,type = "eigen",ecdet = "none",spec = c("longrun"),K=2))##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: maximal eigenvalue statistic (lambda max) , with linear trend
##
## Eigenvalues (lambda):
## [1] 0.23778663 0.16053132 0.11196139 0.02950211
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 3 | 3.53 6.50 8.18 11.65
## r <= 2 | 14.01 12.91 14.90 19.19
## r <= 1 | 20.65 18.90 21.07 25.75
## r = 0 | 32.04 24.78 27.14 32.14
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## logVM.l2 logAM.l2 logSEC.l2 logEXT.l2
## logVM.l2 1.0000000 1.0000000 1.00000000 1.00000000
## logAM.l2 -1.5337623 -1.1270718 0.08946243 0.77995849
## logSEC.l2 -0.4091275 0.3420416 -0.61904255 0.09143807
## logEXT.l2 0.4208877 -0.2642130 -0.79299100 0.28933033
##
## Weights W:
## (This is the loading matrix)
##
## logVM.l2 logAM.l2 logSEC.l2 logEXT.l2
## logVM.d -0.04576371 -0.169638120 -0.05456330 -0.02416623
## logAM.d 0.19703225 0.061873018 -0.04720009 -0.02984947
## logSEC.d 0.57948408 -0.451084842 0.18452091 -0.01718318
## logEXT.d -0.26145793 0.001660205 0.12765052 -0.06088780##Hay evidencia de al menos 1 vector de cointegración al 5%, y posiblemente 2 o 3 al 10%.
#Prueba de JOHANSEN de la traza con constante
summary(ca.jo(datos3,type = "eigen",ecdet = "const",spec = c("longrun"),K=2))##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: maximal eigenvalue statistic (lambda max) , without linear trend and constant in cointegration
##
## Eigenvalues (lambda):
## [1] 2.378559e-01 1.605493e-01 1.120584e-01 3.064810e-02 -1.270190e-18
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 3 | 3.67 7.52 9.24 12.97
## r <= 2 | 14.02 13.75 15.67 20.20
## r <= 1 | 20.65 19.77 22.00 26.81
## r = 0 | 32.05 25.56 28.14 33.24
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## logVM.l2 logAM.l2 logSEC.l2 logEXT.l2 constant
## logVM.l2 1.0000000 1.0000000 1.00000000 1.00000000 1.0000000
## logAM.l2 -1.5340122 -1.1273346 0.08851884 0.76475304 1.2379625
## logSEC.l2 -0.4095895 0.3419297 -0.61826718 0.07754885 0.4730805
## logEXT.l2 0.4209455 -0.2636368 -0.79586595 0.28412769 0.4021763
## constant 0.9224470 1.8922082 3.42289149 -7.90692840 -7.6999659
##
## Weights W:
## (This is the loading matrix)
##
## logVM.l2 logAM.l2 logSEC.l2 logEXT.l2 constant
## logVM.d -0.04552778 -0.169768303 -0.05429868 -0.02453659 1.113079e-16
## logAM.d 0.19713086 0.061864066 -0.04687007 -0.03026915 -6.601442e-17
## logSEC.d 0.57999898 -0.451568236 0.18421321 -0.01690699 1.079698e-16
## logEXT.d -0.26049455 0.000765155 0.12816694 -0.06147254 -1.150006e-17##Hay evidencia de 1 vector de cointegración al 5%, y posiblemente 2 o 3 al 10%.
#Prueba de JOHANSEN de la traza con tendencia
summary(ca.jo(datos3,type = "eigen",ecdet = "trend",spec = c("longrun"),K=2))##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: maximal eigenvalue statistic (lambda max) , with linear trend in cointegration
##
## Eigenvalues (lambda):
## [1] 2.713727e-01 2.201703e-01 1.466493e-01 6.366977e-02 1.665335e-16
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 3 | 7.76 10.49 12.25 16.26
## r <= 2 | 18.71 16.85 18.96 23.65
## r <= 1 | 29.34 23.11 25.54 30.34
## r = 0 | 37.36 29.12 31.46 36.65
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## logVM.l2 logAM.l2 logSEC.l2 logEXT.l2 trend.l2
## logVM.l2 1.000000 1.000000000 1.000000000 1.000000000 1.00000000
## logAM.l2 -13.018481 -1.408017867 -1.291809322 0.460361002 0.20679965
## logSEC.l2 -65.769764 -0.173870988 0.801990168 -0.122778996 0.48830609
## logEXT.l2 197.088088 -0.133256767 0.078772531 -0.079849817 0.15419215
## trend.l2 -1.609663 0.003369843 -0.003784279 -0.006650682 0.03930157
##
## Weights W:
## (This is the loading matrix)
##
## logVM.l2 logAM.l2 logSEC.l2 logEXT.l2 trend.l2
## logVM.d 0.0002490604 -0.17268951 -0.04094915 -0.08554622 -7.741429e-17
## logAM.d 0.0004675289 0.23715128 0.03782472 -0.08490209 -8.010488e-17
## logSEC.d 0.0013239062 0.66119039 -0.37015577 0.02774634 -3.414758e-16
## logEXT.d -0.0024808610 0.09894332 -0.13771592 -0.05503824 5.661135e-16##Hay evidencia de 2 vectores de cointegración al 5%, y posiblemente 3 al 10%.Pronostico
# Estimar VECM
vecm <- ca.jo(datos3, type = "trace", ecdet = "trend", K = 2, spec = "longrun")
vecm_model <- cajorls(vecm, r = 2)# Generar pronósticos
vecm_var <- vec2var(vecm, r = 2)
n_ahead <- 12
vecm_pred <- predict(vecm_var, n.ahead = n_ahead, ci = 0.95)# Transformar a escala original
logVM_pred <- vecm_pred$fcst$logVM[, "fcst"]
logAM_pred <- vecm_pred$fcst$logAM[, "fcst"]
logSEC_pred <- vecm_pred$fcst$logSEC[, "fcst"]
logEXT_pred <- vecm_pred$fcst$logEXT[, "fcst"]
VM_pre <- exp(logVM_pred)
AM_pre <- exp(logAM_pred)
SEC_pre <- exp(logSEC_pred)
EXT_pre <- exp(logEXT_pred)# Crear tabla de predicciones
meses_es <- c("Enero", "Febrero", "Marzo", "Abril", "Mayo", "Junio",
"Julio", "Agosto", "Septiembre", "Octubre", "Noviembre", "Diciembre")
tabla_predicciones <- data.frame(
Mes = meses_es,
Vidrio_manufactura_Proyectada = round(VM_pre, 2),
Aluminio_manufacturas_Proyectada = round(AM_pre, 2),
SECUESTRO_Proyectado = round(SEC_pre, 2),
EXTORSION_Proyectado = round(EXT_pre, 2)
)
print(tabla_predicciones)## Mes Vidrio_manufactura_Proyectada Aluminio_manufacturas_Proyectada
## 1 Enero 26.95 48.06
## 2 Febrero 27.31 47.88
## 3 Marzo 27.36 47.49
## 4 Abril 27.53 47.54
## 5 Mayo 27.60 47.53
## 6 Junio 27.68 47.61
## 7 Julio 27.73 47.68
## 8 Agosto 27.78 47.77
## 9 Septiembre 27.83 47.87
## 10 Octubre 27.87 47.97
## 11 Noviembre 27.91 48.06
## 12 Diciembre 27.96 48.16
## SECUESTRO_Proyectado EXTORSION_Proyectado
## 1 27.86 1299.98
## 2 25.88 1320.95
## 3 25.98 1327.92
## 4 25.69 1332.47
## 5 25.75 1339.51
## 6 25.76 1348.58
## 7 25.84 1359.41
## 8 25.92 1371.25
## 9 26.02 1383.77
## 10 26.11 1396.72
## 11 26.21 1409.97
## 12 26.31 1423.43# Visualizar
par(mfrow = c(2, 2))
plot(1:n_ahead, VM_pre, type = "l", col = "blue", lwd = 2,
main = "Pronóstico: Vidrio_manufactura (VECM)", xlab = "Meses hacia adelante",
ylab = "Vidrio_manufactura", ylim = range(datos$Vidrio_manufactura))
grid()
plot(1:n_ahead, AM_pre, type = "l", col = "blue", lwd = 2,
main = "Pronóstico: Aluminio_manufacturas (VECM)", xlab = "Meses hacia adelante",
ylab = "Aluminio_manufacturas", ylim = range(datos$Aluminio_manufacturas))
grid()
plot(1:n_ahead, SEC_pre, type = "l", col = "blue", lwd = 2,
main = "Pronóstico: SECUESTRO (VECM)", xlab = "Meses hacia adelante",
ylab = "SECUESTRO", ylim = range(datos$SECUESTRO))
grid()
plot(1:n_ahead, EXT_pre, type = "l", col = "blue", lwd = 2,
main = "Pronóstico: EXTORSION (VECM)", xlab = "Meses hacia adelante",
ylab = "EXTORSION", ylim = range(datos$EXTORSION))
grid()
library(readxl)
DATOS_REALES <- read_excel("DATOS_REALES.xlsx")
DATOS_REALES## # A tibble: 3 × 5
## MES Aluminio_manufacturas vidrio_manufactura SECUESTRO EXTORSION
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Enero 63.4 27.1 21 954
## 2 Febrero 49.1 26.5 28 965
## 3 Marzo 48.0 26.2 35 811Comparación
# Crear data frame de pronósticos (primeros 3 meses)
pronosticos <- data.frame(
Mes = c("Enero", "Febrero", "Marzo"),
Vidrio_manufactura_Proyectada = VM_pre[1:3],
Aluminio_manufacturas_Proyectada = AM_pre[1:3],
SECUESTRO_Proyectado = SEC_pre[1:3],
EXTORSION_Proyectado = EXT_pre[1:3]
)
# Datos reales
datos_reales <- data.frame(
Mes = c("Enero", "Febrero", "Marzo"),
Vidrio_manufactura_Real = c(27.130, 26.470, 26.167),
Aluminio_manufacturas_Real = c(63.381, 49.147, 47.985),
SECUESTRO_Real = c(21, 28, 35),
EXTORSION_Real = c(954, 965, 811)
)
# Unir datos
comparacion <- merge(datos_reales, pronosticos, by = "Mes")
print(comparacion)## Mes Vidrio_manufactura_Real Aluminio_manufacturas_Real SECUESTRO_Real
## 1 Enero 27.130 63.381 21
## 2 Febrero 26.470 49.147 28
## 3 Marzo 26.167 47.985 35
## EXTORSION_Real Vidrio_manufactura_Proyectada Aluminio_manufacturas_Proyectada
## 1 954 26.95304 48.05639
## 2 965 27.30556 47.87991
## 3 811 27.36312 47.49251
## SECUESTRO_Proyectado EXTORSION_Proyectado
## 1 27.86236 1299.978
## 2 25.87824 1320.953
## 3 25.97880 1327.920# Calcular errores
calcular_errores <- function(real, pronostico) {
MAE <- mean(abs(real - pronostico))
RMSE <- sqrt(mean((real - pronostico)^2))
return(c(MAE = MAE, RMSE = RMSE))
}
errores_VM <- calcular_errores(comparacion$Vidrio_manufactura_Real, comparacion$Vidrio_manufactura_Proyectada)
errores_AM <- calcular_errores(comparacion$Aluminio_manufacturas_Real, comparacion$Aluminio_manufacturas_Proyectada)
errores_SEC <- calcular_errores(comparacion$SECUESTRO_Real, comparacion$SECUESTRO_Proyectado)
errores_EXT <- calcular_errores(comparacion$EXTORSION_Real, comparacion$EXTORSION_Proyectado)
cat("Errores Vidrio_manufactura:\n", errores_VM, "\n")## Errores Vidrio_manufactura:
## 0.7362134 0.848564cat("Errores Aluminio_manufacturas:\n", errores_AM, "\n")## Errores Aluminio_manufacturas:
## 5.69473 8.882413cat("Errores SECUESTRO:\n", errores_SEC, "\n")## Errores SECUESTRO:
## 6.001773 6.657723cat("Errores EXTORSION:\n", errores_EXT, "\n")## Errores EXTORSION:
## 406.2837 413.7671# Visualizar
par(mfrow = c(2, 2))
plot(1:3, comparacion$Vidrio_manufactura_Real, type = "b", col = "black", lwd = 2, pch = 19,
main = "Vidrio_manufactura: Real vs Pronóstico", xlab = "Mes", ylab = "Vidrio_manufactura",
xaxt = "n", ylim = range(c(comparacion$Vidrio_manufactura_Real, comparacion$Vidrio_manufactura_Proyectada)))
lines(1:3, comparacion$Vidrio_manufactura_Proyectada, type = "b", col = "blue", lwd = 2, pch = 17)
axis(1, at = 1:3, labels = comparacion$Mes)
legend("topright", legend = c("Real", "Pronóstico"), col = c("black", "blue"), pch = c(19, 17), lwd = 2)
grid()
plot(1:3, comparacion$Aluminio_manufacturas_Real, type = "b", col = "black", lwd = 2, pch = 19,
main = "Aluminio_manufacturas: Real vs Pronóstico", xlab = "Mes", ylab = "Aluminio_manufacturas",
xaxt = "n", ylim = range(c(comparacion$Aluminio_manufacturas_Real, comparacion$Aluminio_manufacturas_Proyectada)))
lines(1:3, comparacion$Aluminio_manufacturas_Proyectada, type = "b", col = "blue", lwd = 2, pch = 17)
axis(1, at = 1:3, labels = comparacion$Mes)
legend("topright", legend = c("Real", "Pronóstico"), col = c("black", "blue"), pch = c(19, 17), lwd = 2)
grid()
plot(1:3, comparacion$SECUESTRO_Real, type = "b", col = "black", lwd = 2, pch = 19,
main = "SECUESTRO: Real vs Pronóstico", xlab = "Mes", ylab = "SECUESTRO",
xaxt = "n", ylim = range(c(comparacion$SECUESTRO_Real, comparacion$SECUESTRO_Proyectado)))
lines(1:3, comparacion$SECUESTRO_Proyectado, type = "b", col = "blue", lwd = 2, pch = 17)
axis(1, at = 1:3, labels = comparacion$Mes)
legend("topright", legend = c("Real", "Pronóstico"), col = c("black", "blue"), pch = c(19, 17), lwd = 2)
grid()
plot(1:3, comparacion$EXTORSION_Real, type = "b", col = "black", lwd = 2, pch = 19,
main = "EXTORSION: Real vs Pronóstico", xlab = "Mes", ylab = "EXTORSION",
xaxt = "n", ylim = range(c(comparacion$EXTORSION_Real, comparacion$EXTORSION_Proyectado)))
lines(1:3, comparacion$EXTORSION_Proyectado, type = "b", col = "blue", lwd = 2, pch = 17)
axis(1, at = 1:3, labels = comparacion$Mes)
legend("topright", legend = c("Real", "Pronóstico"), col = c("black", "blue"), pch = c(19, 17), lwd = 2)
grid()
Conclusiones
Las pruebas de Dickey-Fuller (ADF) indicaron que las series en logaritmos (logVM, logAM, logSEC, logEXT) no son estacionarias en niveles, salvo logSEC y logEXT que mostraron estacionariedad con constante y tendencia al 1%, 5% y 10%. Esto justifica el uso de primeras diferencias para el modelo VAR y la exploración de cointegración para VECM.
Las pruebas de Johansen (traza y valor propio máximo) con ecdet = “trend” y K=2 sugieren la existencia de 3 vectores de cointegración al 5% según la prueba de traza, y 2 vectores según la prueba de valor propio máximo. Esto indica relaciones de largo plazo entre las variables, lo que hace adecuado el uso del VECM.
El criterio VARselect sugirió 1 rezago para el VAR y VECM según AIC, HQ, SC y FPE. Sin embargo, el modelo VECM se estimó con K=2 siguiendo las pruebas de cointegración, lo que puede capturar mejor las dinámicas de corto plazo.
Los pronósticos del VECM para Vidrio_manufactura mostraron un MAE de 0.736 y un RMSE de 0.849, indicando una buena precisión, ya que los valores pronosticados (26.95, 27.31, 27.36) están muy cerca de los reales (27.13, 26.47, 26.17). Esto sugiere que el VECM captura bien la dinámica de esta serie.
Los pronósticos del VECM para Vidrio_manufactura y Aluminio_manufacturas están dentro de los rangos históricos (7.7-39.9 y 18.9-73.8, respectivamente), pero los de EXTORSION superan significativamente el rango histórico (88-1354), lo que refuerza la necesidad de revisar la especificación del modelo o incluir variables adicionales para esta serie.