Al igual que en la entrada de Series de Tiempo, este documento es únicamente una guía que me sirve a mí para practicar. No pretendo que funcione como un curso pedagógico para nadie. Si decides que algo te sirve, adelante. Sin embargo, por favor recuerda ir a clases, porque seguramente me equivoqué en algo.

-C







#Vamos a necesitar cargar las siguientes librerías:
library(vars)
## Loading required package: MASS
## 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
## 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)
## Highcharts (www.highcharts.com) is a Highsoft software product which is
## not free for commercial and Governmental use
## 
## Attaching package: 'highcharter'
## The following object is masked from 'package:lmtest':
## 
##     unemployment




A diferencia de los modelos \(AR(p)\) donde las variables endógenas que determinaban el comportamiento futuro de la serie eran retrasos de la serie misma. En los modelos de vectores autoregresivos (\(VAR\)) el comportamiento futuro de la serie será determinado por un vector de variables endógenas. Es decir, el vector de variables endógenas no depende únicamente de los valores pasados de sí mismas, sino de los valores rezagados de otras variables, permitiendo así estimar un sistema de ecuaciones, a diferencia de los modelos ARIMA tradicionales que permiten estimar sólo una ecuación.

Como ejemplo, supón que tenemos tres variables de series de tiempo diferentes, denotadas por \(x_{t,1},x_{t,2},x_{t,3}\). El vector autorregresivo de orden 1 \(VAR(1)\) de la forma reducida de este modelo está dado por \[\begin{align} x_{t,1}&=\alpha_{1}+\beta_{11}x_{t-1,1}+\beta_{12}x_{t-1,2}+\beta_{13}x_{t-1,3}+\varepsilon_{t,1}\\ \\ x_{t,2}&=\alpha_{2}+\beta_{21}x_{t-1,1}+\beta_{22}x_{t-1,2}+\beta_{23}x_{t-1,3}+\varepsilon_{t,2} \tag{1.1} \\ \\ x_{t,3}&=\alpha_{3}+\beta_{31}x_{t-1,1}+\beta_{32}x_{t-1,2}+\beta_{33}x_{t-1,3}+\varepsilon_{t,3} \end{align}\]

Cada variable es una función lineal de los rezagos de un periodo de sí mismas y de las otras dos variables.

Este modelo se puede generalizar a \(n\) variables y expresar matricialmente de la siguiente manera: \[\underbrace{\begin{pmatrix}x_{t,1} \\ x_{t,2} \\ \vdots \\ x_{t,n}\end{pmatrix}}_{X_{t}} = \underbrace{\begin{pmatrix} \alpha_{1} \\ \alpha_{2} \\ \vdots \\ \alpha_{n}\end{pmatrix}}_A + \underbrace{\begin{pmatrix} \beta_{11} & \beta_{12} & \cdots & \beta_{1n}\\ \beta_{21} & \beta_{22} & \cdots & \beta_{2n} \\ \vdots & \vdots & \ddots & \vdots\\ \beta_{n1} & \beta_{n2} & \cdots & \beta_{nn}\end{pmatrix}}_{B} \cdot \underbrace{\begin{pmatrix} x_{t-1,1} \\ x_{t-1,2} \\ \vdots \\ x_{t-1,n} \end{pmatrix}}_{X_{t-1}} + \underbrace{\begin{pmatrix} \varepsilon_{t,1} \\ \varepsilon_{t,2} \\ \vdots \\ \varepsilon_{t,n} \end{pmatrix}}_{E_{t}} \tag{1.2}\] Luego, podemos resumir el modelo general de la ecuación como \(X_{t}=A+BX_{t-1}+E_{t}\) es decir, un sistema de ecuaciones en diferencias de orden 1 con un termino estocastico (\(E_{t}\)).

En general, podemos determinar un \(AR(p)\) rezagado \(p\) con \(n\) variables periodos como \[X_{t}=A+\sum_{i=1}^{p}B_{i}X_{t-i}+E_{t}\] Donde cada \(B_{i}\) es la matriz de \(n\times n\) coeficientes asociadas al periodo \(t-i\).

Construcción del Modelo

Debemos seguir el procedimiento iterativo de Box-Jenkins:

  1. \(\textbf{Identificación}:\) En lugar de identificar un solo modelo ARIMA, debes identificar las relaciones entre las series de tiempo y determinar cuántos rezagos (lags) son necesarios para modelar adecuadamente las dependencias temporales. Esto implica analizar las funciones de autocorrelación y autocorrelación parcial de cada serie de tiempo, así como las funciones de autocorrelación cruzada entre las series.

El primer paso en el procedimiento iterativo de Box-Jenkins para el análisis de series de tiempo implica la identificación de un modelo apropiado. Esto incluye la determinación de los rezagos (lags) necesarios y la comprensión de la estructura de la serie de tiempo.

Cargar y preparar los datos:

data <- as.data.frame(Canada)                        #Esta es una base de datos de la librería vars
colnames(data) <- c("Employment", "Productivity",
                    "Real Wage", "Unemployment Rate")
data <- ts(data, start = c(1980,1),frequency = 4)

#Gráfico de las Series
par(mfrow = c(2, 2))
ts.plot(data[,1], 
        main = colnames(data)[1],
        ylab = "Valor")
ts.plot(data[,2], 
        main = colnames(data)[2],
        ylab = "Valor")
ts.plot(data[,3], 
        main = colnames(data)[3],
        ylab = "Valor")
ts.plot(data[,4], 
        main = colnames(data)[4],
        ylab = "Valor")

Pruebas de estacionariedad:

Al igual que en los modelos univariados para utilizar la metodología \(VAR\) para el modelamiento y posterior pronóstico de series de tiempo se debe verificar la estacionaridad de la serie. Esta prueba se puede hacer a través de un test de Dickey-Fuller (adf.test()) o con la prueba de Phillips-Perron (pp.test()) serie a serie.

adf.test(data[,1])
## 
##  Augmented Dickey-Fuller Test
## 
## data:  data[, 1]
## Dickey-Fuller = -2.148, Lag order = 4, p-value = 0.5152
## alternative hypothesis: stationary
pp.test(data[,1])
## 
##  Phillips-Perron Unit Root Test
## 
## data:  data[, 1]
## Dickey-Fuller Z(alpha) = -5.8701, Truncation lag parameter = 3, p-value
## = 0.7735
## alternative hypothesis: stationary

Ambas pruebas arrojan un \(p-value\) muy elevado, por lo tanto podemos decir que la serie de empleo es estacionaria, por lo que habrá que diferenciar la serie en el siguiente paso.

adf.test(data[,2])
## 
##  Augmented Dickey-Fuller Test
## 
## data:  data[, 2]
## Dickey-Fuller = -2.8725, Lag order = 4, p-value = 0.218
## alternative hypothesis: stationary
pp.test(data[,2])
## 
##  Phillips-Perron Unit Root Test
## 
## data:  data[, 2]
## Dickey-Fuller Z(alpha) = -7.4247, Truncation lag parameter = 3, p-value
## = 0.6816
## alternative hypothesis: stationary
adf.test(data[,3])
## 
##  Augmented Dickey-Fuller Test
## 
## data:  data[, 3]
## Dickey-Fuller = -2.0558, Lag order = 4, p-value = 0.553
## alternative hypothesis: stationary
pp.test(data[,3])
## 
##  Phillips-Perron Unit Root Test
## 
## data:  data[, 3]
## Dickey-Fuller Z(alpha) = -4.6296, Truncation lag parameter = 3, p-value
## = 0.8468
## alternative hypothesis: stationary
adf.test(data[,4])
## 
##  Augmented Dickey-Fuller Test
## 
## data:  data[, 4]
## Dickey-Fuller = -2.5988, Lag order = 4, p-value = 0.3303
## alternative hypothesis: stationary
pp.test(data[,4])
## 
##  Phillips-Perron Unit Root Test
## 
## data:  data[, 4]
## Dickey-Fuller Z(alpha) = -7.1861, Truncation lag parameter = 3, p-value
## = 0.6957
## alternative hypothesis: stationary

Ya que ninguna serie paso el test de estacionalidad tenemos que diferenciar las series mediante la funcion diff()

data.diff <- diff(data[,c(1:4)])

#Empleo
adf.test(data.diff[,1])
## 
##  Augmented Dickey-Fuller Test
## 
## data:  data.diff[, 1]
## Dickey-Fuller = -3.2668, Lag order = 4, p-value = 0.08274
## alternative hypothesis: stationary
pp.test(data.diff[,1])
## 
##  Phillips-Perron Unit Root Test
## 
## data:  data.diff[, 1]
## Dickey-Fuller Z(alpha) = -26.227, Truncation lag parameter = 3, p-value
## = 0.01233
## alternative hypothesis: stationary
#Productividad
adf.test(data.diff[,2])
## 
##  Augmented Dickey-Fuller Test
## 
## data:  data.diff[, 2]
## Dickey-Fuller = -3.2361, Lag order = 4, p-value = 0.08775
## alternative hypothesis: stationary
pp.test(data.diff[,2])
## Warning in pp.test(data.diff[, 2]): p-value smaller than printed p-value
## 
##  Phillips-Perron Unit Root Test
## 
## data:  data.diff[, 2]
## Dickey-Fuller Z(alpha) = -60.867, Truncation lag parameter = 3, p-value
## = 0.01
## alternative hypothesis: stationary
#Salario Real
adf.test(data.diff[,3])
## 
##  Augmented Dickey-Fuller Test
## 
## data:  data.diff[, 3]
## Dickey-Fuller = -3.2352, Lag order = 4, p-value = 0.08789
## alternative hypothesis: stationary
pp.test(data.diff[,3])
## Warning in pp.test(data.diff[, 3]): p-value smaller than printed p-value
## 
##  Phillips-Perron Unit Root Test
## 
## data:  data.diff[, 3]
## Dickey-Fuller Z(alpha) = -61.576, Truncation lag parameter = 3, p-value
## = 0.01
## alternative hypothesis: stationary
#Tasa de Desempleo
adf.test(data.diff[,4])
## 
##  Augmented Dickey-Fuller Test
## 
## data:  data.diff[, 4]
## Dickey-Fuller = -3.7325, Lag order = 4, p-value = 0.02698
## alternative hypothesis: stationary
pp.test(data.diff[,4])
## Warning in pp.test(data.diff[, 4]): p-value smaller than printed p-value
## 
##  Phillips-Perron Unit Root Test
## 
## data:  data.diff[, 4]
## Dickey-Fuller Z(alpha) = -38.946, Truncation lag parameter = 3, p-value
## = 0.01
## alternative hypothesis: stationary

Gráfica de las series diferenciadas:

par(mfrow = c(2, 2))
ts.plot(data.diff[,1], 
        main = colnames(data.diff)[1],
        ylab = "Valor")
ts.plot(data.diff[,2], 
        main = colnames(data.diff)[2],
        ylab = "Valor")
ts.plot(data.diff[,3], 
        main = colnames(data.diff)[3],
        ylab = "Valor")
ts.plot(data.diff[,4], 
        main = colnames(data.diff)[4],
        ylab = "Valor")

Estimación:

Una vez que tenemos las series estacionarias procedemos a hacer las estimaciones correspondientes.

## Debemos conocer el numero de rezagos apropiados para nuestro modelo:
VARselect(data.diff, lag.max = 8,
          type = "const")         #Esta función nos va a permitir conocer el numero de rezagos apropiados para simular nuesto modelo.
## $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
## Generar el modelo mediante la función VAR():
modelo_var<-VAR(data, p=2)        #El valor p=1 es el número de rezagos lo que modelará un VAR(1)
summary(modelo_var)
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: Employment, Productivity, Real.Wage, Unemployment.Rate 
## 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 = data, p = 2)
## 
## 
## Estimation results for equation Employment: 
## =========================================== 
## Employment = Employment.l1 + Productivity.l1 + Real.Wage.l1 + Unemployment.Rate.l1 + Employment.l2 + Productivity.l2 + Real.Wage.l2 + Unemployment.Rate.l2 + const 
## 
##                        Estimate Std. Error t value Pr(>|t|)    
## Employment.l1         1.638e+00  1.500e-01  10.918  < 2e-16 ***
## Productivity.l1       1.673e-01  6.114e-02   2.736  0.00780 ** 
## Real.Wage.l1         -6.312e-02  5.524e-02  -1.143  0.25692    
## Unemployment.Rate.l1  2.656e-01  2.028e-01   1.310  0.19444    
## Employment.l2        -4.971e-01  1.595e-01  -3.116  0.00262 ** 
## Productivity.l2      -1.017e-01  6.607e-02  -1.539  0.12824    
## Real.Wage.l2          3.844e-03  5.552e-02   0.069  0.94499    
## Unemployment.Rate.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 Productivity: 
## ============================================= 
## Productivity = Employment.l1 + Productivity.l1 + Real.Wage.l1 + Unemployment.Rate.l1 + Employment.l2 + Productivity.l2 + Real.Wage.l2 + Unemployment.Rate.l2 + const 
## 
##                        Estimate Std. Error t value Pr(>|t|)    
## Employment.l1          -0.17277    0.26977  -0.640  0.52390    
## Productivity.l1         1.15043    0.10995  10.464 3.57e-16 ***
## Real.Wage.l1            0.05130    0.09934   0.516  0.60710    
## Unemployment.Rate.l1   -0.47850    0.36470  -1.312  0.19362    
## Employment.l2           0.38526    0.28688   1.343  0.18346    
## Productivity.l2        -0.17241    0.11881  -1.451  0.15104    
## Real.Wage.l2           -0.11885    0.09985  -1.190  0.23778    
## Unemployment.Rate.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 Real.Wage: 
## ========================================== 
## Real.Wage = Employment.l1 + Productivity.l1 + Real.Wage.l1 + Unemployment.Rate.l1 + Employment.l2 + Productivity.l2 + Real.Wage.l2 + Unemployment.Rate.l2 + const 
## 
##                        Estimate Std. Error t value Pr(>|t|)    
## Employment.l1         -0.268833   0.322619  -0.833    0.407    
## Productivity.l1       -0.081065   0.131487  -0.617    0.539    
## Real.Wage.l1           0.895478   0.118800   7.538 1.04e-10 ***
## Unemployment.Rate.l1   0.012130   0.436149   0.028    0.978    
## Employment.l2          0.367849   0.343087   1.072    0.287    
## Productivity.l2       -0.005181   0.142093  -0.036    0.971    
## Real.Wage.l2           0.052677   0.119410   0.441    0.660    
## Unemployment.Rate.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 Unemployment.Rate: 
## ================================================== 
## Unemployment.Rate = Employment.l1 + Productivity.l1 + Real.Wage.l1 + Unemployment.Rate.l1 + Employment.l2 + Productivity.l2 + Real.Wage.l2 + Unemployment.Rate.l2 + const 
## 
##                       Estimate Std. Error t value Pr(>|t|)    
## Employment.l1         -0.58076    0.11563  -5.023 3.49e-06 ***
## Productivity.l1       -0.07812    0.04713  -1.658 0.101682    
## Real.Wage.l1           0.01866    0.04258   0.438 0.662463    
## Unemployment.Rate.l1   0.61893    0.15632   3.959 0.000173 ***
## Employment.l2          0.40982    0.12296   3.333 0.001352 ** 
## Productivity.l2        0.05212    0.05093   1.023 0.309513    
## Real.Wage.l2           0.04180    0.04280   0.977 0.331928    
## Unemployment.Rate.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:
##                   Employment Productivity Real.Wage Unemployment.Rate
## Employment          0.131635    -0.007469  -0.04210          -0.06909
## Productivity       -0.007469     0.425711   0.06461           0.01392
## Real.Wage          -0.042099     0.064613   0.60886           0.03422
## Unemployment.Rate  -0.069087     0.013923   0.03422           0.07821
## 
## Correlation matrix of residuals:
##                   Employment Productivity Real.Wage Unemployment.Rate
## Employment           1.00000     -0.03155   -0.1487           -0.6809
## Productivity        -0.03155      1.00000    0.1269            0.0763
## Real.Wage           -0.14870      0.12691    1.0000            0.1568
## Unemployment.Rate   -0.68090      0.07630    0.1568            1.0000

La elección del \(VAR(1)\) fue debido a que los valores \(p\) asociados a las regresiones eran estadisticamente significativos cuando utilizabamos rezagos más elevados.

Gráficamos el modelo:

plot(modelo_var)

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

Al ser todas menores a 1 este modelo cumple el supuesto de estacionieradad.

Comprobación de la autocorrelación en los errores:

La comprobación de la autocorrelación en los errores es una parte importante del análisis de modelos de series temporales y regresión. La autocorrelación se refiere a la correlación de un valor en una serie temporal o una secuencia de datos con sus valores anteriores o posteriores en la misma serie. En el contexto de los modelos de regresión, la autocorrelación en los errores puede indicar que el modelo no está capturando adecuadamente la estructura temporal de los datos o que los errores no son independientes y, por lo tanto, no se cumplen los supuestos básicos del modelo

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

La función serial.test() se utiliza para realizar prueba de Ljung-Box para evaluar la autocorrelación en los residuos del modelo. Estas pruebas evalúan si hay correlación en los residuos en diferentes retrasos (lags).

El resultado de la prueba será un valor p (p-value) que indica la probabilidad de observar la autocorrelación en los residuos si los residuos fueran realmente independientes y no autocorrelacionados. En general, si el valor p es pequeño (por ejemplo, menor que 0.05), se considera evidencia de autocorrelación significativa en los residuos, lo que sugiere que el modelo puede no ser adecuado.

Por lo tanto, en este ejemplo, no hay autocorrelacion de los errores

Prueba de nomalidad:

normality.test(modelo_var,
               multivariate.only = F)
## $Employment
## 
##  JB-Test (univariate)
## 
## data:  Residual of Employment equation
## Chi-squared = 0.15347, df = 2, p-value = 0.9261
## 
## 
## $Productivity
## 
##  JB-Test (univariate)
## 
## data:  Residual of Productivity equation
## Chi-squared = 4.2651, df = 2, p-value = 0.1185
## 
## 
## $Real.Wage
## 
##  JB-Test (univariate)
## 
## data:  Residual of Real.Wage equation
## Chi-squared = 0.33482, df = 2, p-value = 0.8459
## 
## 
## $Unemployment.Rate
## 
##  JB-Test (univariate)
## 
## data:  Residual of Unemployment.Rate 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
  • Si p.value es menor que el nivel de significancia (0.05), se rechaza la hipótesis nula y se concluye que los residuos no siguen una distribución normal.
  • Si p.value es mayor o igual al nivel de significancia, no hay suficiente evidencia para afirmar que los residuos no siguen una distribución normal.
residuos <- residuals(modelo_var)

par(mfrow = c(2, 1))

hist(residuos, main = "Histograma de Residuos", xlab = "Valor Residual", prob = TRUE, col = "lightblue")


mu <- mean(residuos)
sigma <- sd(residuos)
x <- seq(min(residuos), max(residuos), length = 100)
y <- dnorm(x, mean = mu, sd = sigma)
lines(x, y, col = "red", lwd = 2)


legend("topright", legend = "Densidad Normal", col = "red", lty = 1, lwd = 2)

# Crear un gráfico QQ para evaluar la normalidad
qqnorm(residuos)
qqline(residuos)

Pronostico:

pronosticos <- predict(modelo_var,
                       n.ahead = 2)                    #El pronostico del modelo_var se realizará para                                                        #dos trimestres adelante para cada una de las series.

## Pronosticos para las series:

pronostico_e <- pronosticos$fcst$Employment            #Pronostico de la serie de empleo

pronostico_p <- pronosticos$fcst$Productivity          #Pronostico de la serie de productividad

pronostico_rw <- pronosticos$fcst$Real.Wage.           #Pronostico de la serie de salario real

pronostico_ur <- pronosticos$fcst$Unemployment.Rate.   #Pronostico de la serie de tasa de desempleo

## Intervalos de condianza de cada pronostico:

# Empleo
intervalos_inf_e <- pronosticos$fcst$Employment[,2]
intervalos_sup_e <- pronosticos$fcst$Employment[,3]

# Productividad
intervalos_inf_p <- pronosticos$fcst$Productivity[,2]
intervalos_sup_p <- pronosticos$fcst$Productivity[,3]

# Salario Real
intervalos_inf_rw <- pronosticos$fcst$Real.Wage[,2]
intervalos_sup_rw <- pronosticos$fcst$Real.Wage[,3]

# Tasa de desempleo
intervalos_inf_ur <- pronosticos$fcst$Unemployment.Rate[,2]
intervalos_sup_ur <- pronosticos$fcst$Unemployment.Rate[,3]







Cuiudate bb