Aplicación 3.6: Autocorrelación temporal y regresiones dinámicas

Caso univariante (especificación ARDL): Tipos de interés en España

En la primera aplicación se estimará un modelo de regresión uniecuacional para explicar los tipos de interés de largo plazo vigentes en la economía española en la década de los años 1980, se analizará la correlación serial de los residuos de dicho modelo, y se propondrán soluciones al problema detectado de autocorrelación en los errores.

La regresión de referencia utilizada es la siguiente:

\[RL_{t} = \beta_1 + \beta_2 R3M_{t} + \beta_3 RD_{t} + e_{t}\]

donde \(RL\) y \(R3M\) son los tipos de interés de largo y corto (3 meses) plazo, respectivamente, y \(RD\) es el tipo de interés de remuneración de los depósitos bancarios. Dicha especificación se corresponde con el denominado modelo de Klein-Monti.

Código R

# Lectura de librerías
library(tidyverse)
library(dynlm)
library(car)
library(lmtest)
library(sandwich)
library(orcutt)
# Lectura de datos
TIPOS_INT_ESP <- read_delim("data/TIPOS_INT_ESP.csv", delim = ";")
str(TIPOS_INT_ESP)

spc_tbl_ [35 × 4] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
 $ date: Date[1:35], format: "1982-01-01" "1982-04-01" ...
 $ R3M : num [1:35] 16.1 17.5 17.1 18.5 17.9 ...
 $ RD  : num [1:35] 8.12 8.29 8.36 8.49 8.54 8.7 8.71 8.79 9.04 9.12 ...
 $ RL  : num [1:35] 17.5 17.6 17.7 17.7 17.7 ...
 - attr(*, "spec")=
  .. cols(
  ..   date = col_date(format = ""),
  ..   R3M = col_double(),
  ..   RD = col_double(),
  ..   RL = col_double()
  .. )
 - attr(*, "problems")=<externalptr> 

# Gráfica de las series temporales
ggplot(TIPOS_INT_ESP,aes(date,R3M)) +
  geom_line(color="blue") +
  geom_line(aes(date,RD),color="orange") +
  geom_line(aes(date,RL),color="green") +
  xlab("") + ylab("Tipos de interés")

# Asignación del formato trimestral
TIPOS_ESP_ts <- ts(TIPOS_INT_ESP[,2:4], start=c(1982,1), end = c(1990,3), 
                   frequency = 4)
plot(TIPOS_ESP_ts)

# Modelo de regresión estático
# Función lm() estándar:
# summary(lm_KM <- lm(RL ~ R3M + RD, data=TIPOS_ESP_ts))
# Librería dynlm (especificación ARDL(0,0,0))
summary(dynlm_KM_0 <- dynlm(RL ~ R3M + RD, data=TIPOS_ESP_ts)) 

Time series regression with "ts" data:
Start = 1982(1), End = 1990(3)

Call:
dynlm(formula = RL ~ R3M + RD, data = TIPOS_ESP_ts)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.75780 -0.16274  0.07166  0.25488  0.61088 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  7.02728    0.48653  14.444 1.43e-15 ***
R3M          0.15918    0.02301   6.918 7.84e-08 ***
RD           0.92600    0.07777  11.907 2.69e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3567 on 32 degrees of freedom
Multiple R-squared:  0.9322,    Adjusted R-squared:  0.9279 
F-statistic: 219.9 on 2 and 32 DF,  p-value: < 2.2e-16

# Contrastes de correlación en los errores (autocorrelación)
# Residuos del modelo
resid <-residuals(dynlm_KM_0)
plot(resid)
abline(h=0, lty=2)

# Correlograma de los residuos
corr <- acf(resid)

corr$acf[2:10]

[1] 0.661182301 0.373478408 0.196080083 0.059991302 0.055270320 0.096234644
[7] 0.258730586 0.172078454 0.005038669

# Test de Durbin-Watson
dwtest(dynlm_KM_0, alternative = "two.sided")

    Durbin-Watson test

data:  dynlm_KM_0
DW = 0.49682, p-value = 9.951e-09
alternative hypothesis: true autocorrelation is not 0

dwtest(dynlm_KM_0, alternative = "greater")

    Durbin-Watson test

data:  dynlm_KM_0
DW = 0.49682, p-value = 4.976e-09
alternative hypothesis: true autocorrelation is greater than 0

# Test de Breusch-Godfrey
bgtest(dynlm_KM_0, order=1, type="Chisq", fill=0)

    Breusch-Godfrey test for serial correlation of order up to 1

data:  dynlm_KM_0
LM test = 18.689, df = 1, p-value = 1.539e-05

# Corrección de la autocorrelación
# MCO corregidos: errores estándar robustos, HAC (Newey-West)
summary(dynlm_KM_0 <- dynlm(RL ~ R3M + RD, data=TIPOS_ESP_ts), 
        vcov.=vcovHAC(dynlm_KM_0))

Time series regression with "ts" data:
Start = 1982(1), End = 1990(3)

Call:
dynlm(formula = RL ~ R3M + RD, data = TIPOS_ESP_ts)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.75780 -0.16274  0.07166  0.25488  0.61088 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  7.02728    0.66475  10.571 5.72e-12 ***
R3M          0.15918    0.03572   4.457 9.56e-05 ***
RD           0.92600    0.11776   7.863 5.69e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3567 on 32 degrees of freedom
Multiple R-squared:  0.9322,    Adjusted R-squared:  0.9279 
F-statistic: 146.4 on 2 and 32 DF,  p-value: < 2.2e-16

# Mínimos cuadrados generalizados (MCG): errores AR(1)
cochrane.orcutt(dynlm_KM_0)

Cochrane-orcutt estimation for first order autocorrelation 
 
Call:
dynlm(formula = RL ~ R3M + RD, data = TIPOS_ESP_ts)

 number of interaction: 11
 rho 0.773176

Durbin-Watson statistic 
(original):    0.49682 , p-value: 4.976e-09
(transformed): 1.69696 , p-value: 1.308e-01
 
 coefficients: 
(Intercept)         R3M          RD 
   7.013488    0.157823    0.913996 

# Dinámica en las variables: modelo ARDL(1,1,1) 
summary(dynlm_KM_1 <- dynlm(RL ~ L(RL, 1:1) + L(R3M, 0:1) + L(RD, 0:1), 
                            data=TIPOS_ESP_ts))

Time series regression with "ts" data:
Start = 1982(2), End = 1990(3)

Call:
dynlm(formula = RL ~ L(RL, 1:1) + L(R3M, 0:1) + L(RD, 0:1), data = TIPOS_ESP_ts)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.40953 -0.14242 -0.01385  0.16811  0.38224 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)   1.95931    1.02773   1.906  0.06690 .  
L(RL, 1:1)    0.71922    0.13830   5.200 1.60e-05 ***
L(R3M, 0:1)0  0.16624    0.02491   6.673 3.05e-07 ***
L(R3M, 0:1)1 -0.07777    0.03656  -2.127  0.04237 *  
L(RD, 0:1)0   0.65959    0.19117   3.450  0.00179 ** 
L(RD, 0:1)1  -0.48983    0.24726  -1.981  0.05748 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2333 on 28 degrees of freedom
Multiple R-squared:  0.9743,    Adjusted R-squared:  0.9697 
F-statistic:   212 on 5 and 28 DF,  p-value: < 2.2e-16

# Comparación de modelos
compareCoefs(dynlm_KM_0,dynlm_KM_1)

Calls:
1: dynlm(formula = RL ~ R3M + RD, data = TIPOS_ESP_ts)
2: dynlm(formula = RL ~ L(RL, 1:1) + L(R3M, 0:1) + L(RD, 0:1), data = 
  TIPOS_ESP_ts)

             Model 1 Model 2
(Intercept)    7.027   1.959
SE             0.487   1.028
                            
R3M            0.159        
SE             0.023        
                            
RD            0.9260        
SE            0.0778        
                            
L(RL, 1:1)             0.719
SE                     0.138
                            
L(R3M, 0:1)0          0.1662
SE                    0.0249
                            
L(R3M, 0:1)1         -0.0778
SE                    0.0366
                            
L(RD, 0:1)0            0.660
SE                     0.191
                            
L(RD, 0:1)1           -0.490
SE                     0.247
                            

# Estimación de los efectos parciales asociados
# Efectos a corto y largo plazo
library(nlWaldTest)
# Modelo estático (cp=lp)
nlConfint(dynlm_KM_0, c("b[2]","b[3]"))

         value     2.5 %    97.5 %
b[2] 0.1591826 0.1140841 0.2042812
b[3] 0.9260003 0.7735733 1.0784272

nlWaldtest(dynlm_KM_0, "b[2]")

    Wald Chi-square test of a restriction on model parameters

data:  dynlm_KM_0
Chisq = 47.859, df = 1, p-value = 4.58e-12

nlWaldtest(dynlm_KM_0, "b[3]")

    Wald Chi-square test of a restriction on model parameters

data:  dynlm_KM_0
Chisq = 141.77, df = 1, p-value < 2.2e-16

# Modelo dinámico
# Corto plazo
nlConfint(dynlm_KM_1, c("b[3]","b[5]"))

         value     2.5 %    97.5 %
b[3] 0.1662367 0.1174099 0.2150636
b[5] 0.6595859 0.2848907 1.0342812

nlWaldtest(dynlm_KM_1, "b[3]")

    Wald Chi-square test of a restriction on model parameters

data:  dynlm_KM_1
Chisq = 44.528, df = 1, p-value = 2.508e-11

nlWaldtest(dynlm_KM_1, "b[5]")

    Wald Chi-square test of a restriction on model parameters

data:  dynlm_KM_1
Chisq = 11.904, df = 1, p-value = 0.0005602

# Largo plazo
nlConfint(dynlm_KM_1, c("(b[3]+b[4])/(1-b[2])","(b[5]+b[6])/(1-b[2])"))

                         value     2.5 %    97.5 %
(b[3]+b[4])/(1-b[2]) 0.3150717 0.1198835 0.5102598
(b[5]+b[6])/(1-b[2]) 0.6045916 0.1048464 1.1043367

nlWaldtest(dynlm_KM_1, "(b[3]+b[4])/(1-b[2])")

    Wald Chi-square test of a restriction on model parameters

data:  dynlm_KM_1
Chisq = 10.009, df = 1, p-value = 0.001557

nlWaldtest(dynlm_KM_1, "(b[5]+b[6])/(1-b[2])")

    Wald Chi-square test of a restriction on model parameters

data:  dynlm_KM_1
Chisq = 5.6224, df = 1, p-value = 0.01773

Código Python

# Lectura de librerías
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import statsmodels.api as sm
import statsmodels.formula.api as smf
import statsmodels.stats.api as sms
from statsmodels.graphics import tsaplots
import statsmodels.stats.diagnostic as smsdiag
from statsmodels.compat import lzip
# Lectura de datos y asignación del índice temporal
TIPOS_INT_ESP = pd.read_csv("data/TIPOS_INT_ESP.csv", sep=";", 
parse_dates=['date'], index_col='date')
TIPOS_INT_ESP.info()

<class 'pandas.core.frame.DataFrame'>
DatetimeIndex: 35 entries, 1982-01-01 to 1990-07-01
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   R3M     35 non-null     float64
 1   RD      35 non-null     float64
 2   RL      35 non-null     float64
dtypes: float64(3)
memory usage: 1.1 KB

# Gráfica de las series temporales
plt.figure(1)
fig, ax = plt.subplots()
TIPOS_INT_ESP.plot(ax=ax)
plt.legend(['R3M','RD','RL']); plt.xlabel(''); plt.ylabel('Tipos de interés')
plt.show()

# Asignación del formato trimestral
TIPOS_ESP_ts = pd.read_csv("data/TIPOS_INT_ESP.csv", sep=";").iloc[:,1:]
date = pd.date_range(start = '1982', periods = len(TIPOS_ESP_ts.index), 
freq = 'QS')
TIPOS_ESP_ts.index = date
TIPOS_ESP_ts.info()

<class 'pandas.core.frame.DataFrame'>
DatetimeIndex: 35 entries, 1982-01-01 to 1990-07-01
Freq: QS-JAN
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   R3M     35 non-null     float64
 1   RD      35 non-null     float64
 2   RL      35 non-null     float64
dtypes: float64(3)
memory usage: 1.1 KB

TIPOS_ESP_ts.index

DatetimeIndex(['1982-01-01', '1982-04-01', '1982-07-01', '1982-10-01',
               '1983-01-01', '1983-04-01', '1983-07-01', '1983-10-01',
               '1984-01-01', '1984-04-01', '1984-07-01', '1984-10-01',
               '1985-01-01', '1985-04-01', '1985-07-01', '1985-10-01',
               '1986-01-01', '1986-04-01', '1986-07-01', '1986-10-01',
               '1987-01-01', '1987-04-01', '1987-07-01', '1987-10-01',
               '1988-01-01', '1988-04-01', '1988-07-01', '1988-10-01',
               '1989-01-01', '1989-04-01', '1989-07-01', '1989-10-01',
               '1990-01-01', '1990-04-01', '1990-07-01'],
              dtype='datetime64[ns]', freq='QS-JAN')

plt.figure(2)
fig, ax = plt.subplots(3, 1, sharex = True)
ax[0].plot(TIPOS_ESP_ts.R3M)
ax[0].set_ylabel('R3M')
ax[1].plot(TIPOS_ESP_ts.RD)
ax[1].set_ylabel('RD')
ax[2].plot(TIPOS_ESP_ts.RL)
ax[2].set_ylabel('RL')
plt.show()

# Modelo de regresión estático
model = smf.ols(formula = "RL ~ R3M + RD", data = TIPOS_ESP_ts)
lm_KM = model.fit()
print(lm_KM.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                     RL   R-squared:                       0.932
Model:                            OLS   Adj. R-squared:                  0.928
Method:                 Least Squares   F-statistic:                     219.9
Date:                Sun, 09 Feb 2025   Prob (F-statistic):           2.00e-19
Time:                        14:21:36   Log-Likelihood:                -12.013
No. Observations:                  35   AIC:                             30.03
Df Residuals:                      32   BIC:                             34.69
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept      7.0273      0.487     14.444      0.000       6.036       8.018
R3M            0.1592      0.023      6.918      0.000       0.112       0.206
RD             0.9260      0.078     11.907      0.000       0.768       1.084
==============================================================================
Omnibus:                        3.138   Durbin-Watson:                   0.497
Prob(Omnibus):                  0.208   Jarque-Bera (JB):                2.855
Skew:                          -0.671   Prob(JB):                        0.240
Kurtosis:                       2.602   Cond. No.                         144.
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

# Contrastes de correlación en los errores (autocorrelación)
# Errores estimados (residuos)
residuos = lm_KM.resid
plt.figure(3)
plt.plot(residuos, label='Residuos')
plt.xlabel('Date')
plt.legend()
plt.show()

# Correlograma de los residuos
r = sm.tsa.stattools.acf(residuos, nlags=9, fft=True)
plt.figure(4)
fig, ax = plt.subplots()
tsaplots.plot_acf(residuos, lags=15, alpha=0.05, zero=False, 
vlines_kwargs={"colors":'black'}, color='black', title='', ax=ax)
plt.show()

print(r)

[1.         0.6611823  0.37347841 0.19608008 0.0599913  0.05527032
 0.09623464 0.25873059 0.17207845 0.00503867]

# Test de Durbin-Watson
sms.durbin_watson(residuos).round(5)

0.49682

# Test de Breusch-Godfrey
name = ['LM statistic', 'Chi^2 p-val', 'F statistic', 'F p-val']
BG_test = smsdiag.acorr_breusch_godfrey(lm_KM,1)
lzip(name, BG_test)

[('LM statistic', 18.689051294910442), ('Chi^2 p-val', 1.538634377813448e-05), ('F statistic', 35.51973589135522), ('F p-val', 1.3756949593855872e-06)]

# Corrección de la autocorrelación
# MCO corregidos: errores estándar robustos
# Corrección de la matriz de covarianzas: estimador de Newey-West 
# (si maxlags=0 -> HAC=HC1)
lm_KM_HAC = smf.ols(formula = "RL ~ R3M + RD", data = TIPOS_ESP_ts).fit(cov_type='HAC', cov_kwds={'maxlags':6,'use_correction':True})
print(lm_KM_HAC.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                     RL   R-squared:                       0.932
Model:                            OLS   Adj. R-squared:                  0.928
Method:                 Least Squares   F-statistic:                     103.2
Date:                Sun, 09 Feb 2025   Prob (F-statistic):           1.11e-14
Time:                        14:21:37   Log-Likelihood:                -12.013
No. Observations:                  35   AIC:                             30.03
Df Residuals:                      32   BIC:                             34.69
Df Model:                           2                                         
Covariance Type:                  HAC                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept      7.0273      0.859      8.183      0.000       5.344       8.710
R3M            0.1592      0.033      4.778      0.000       0.094       0.224
RD             0.9260      0.142      6.509      0.000       0.647       1.205
==============================================================================
Omnibus:                        3.138   Durbin-Watson:                   0.497
Prob(Omnibus):                  0.208   Jarque-Bera (JB):                2.855
Skew:                          -0.671   Prob(JB):                        0.240
Kurtosis:                       2.602   Cond. No.                         144.
==============================================================================

Notes:
[1] Standard Errors are heteroscedasticity and autocorrelation robust (HAC) using 6 lags and with small sample correction

# Dinámica en las variables: modelo ARDL(1,1,1) 
model = smf.ols(formula = "RL ~ RL.shift(1) + R3M + R3M.shift(1) + RD + RD.shift(1)", 
data = TIPOS_ESP_ts)
lm_KM_dyn = model.fit()
print(lm_KM_dyn.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                     RL   R-squared:                       0.974
Model:                            OLS   Adj. R-squared:                  0.970
Method:                 Least Squares   F-statistic:                     212.0
Date:                Sun, 09 Feb 2025   Prob (F-statistic):           2.41e-21
Time:                        14:21:37   Log-Likelihood:                 4.5383
No. Observations:                  34   AIC:                             2.923
Df Residuals:                      28   BIC:                             12.08
Df Model:                           5                                         
Covariance Type:            nonrobust                                         
================================================================================
                   coef    std err          t      P>|t|      [0.025      0.975]
--------------------------------------------------------------------------------
Intercept        1.9593      1.028      1.906      0.067      -0.146       4.065
RL.shift(1)      0.7192      0.138      5.200      0.000       0.436       1.003
R3M              0.1662      0.025      6.673      0.000       0.115       0.217
R3M.shift(1)    -0.0778      0.037     -2.127      0.042      -0.153      -0.003
RD               0.6596      0.191      3.450      0.002       0.268       1.051
RD.shift(1)     -0.4898      0.247     -1.981      0.057      -0.996       0.017
==============================================================================
Omnibus:                        2.342   Durbin-Watson:                   1.600
Prob(Omnibus):                  0.310   Jarque-Bera (JB):                1.374
Skew:                          -0.177   Prob(JB):                        0.503
Kurtosis:                       2.081   Cond. No.                         787.
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

b = lm_KM_dyn.params[1:]
# Multiplicadores de corto plazo
b[1].round(3) ; b[3].round(3)

0.166
0.66

# Multiplicadores de largo plazo
b_lr_R3M = (b[1]+b[2])/(1-b[0])
b_lr_R3M.round(3)

0.315

b_lr_RD = (b[3]+b[4])/(1-b[0])
b_lr_RD.round(3)

0.605

# Contrastes
H_0_b1 = '(R3M = 0)'
W_test = lm_KM_dyn.wald_test(H_0_b1, use_f = False)
W_test

<class 'statsmodels.stats.contrast.ContrastResults'>
<Wald test (chi2): statistic=[[44.52791871]], p-value=2.507557561604631e-11, df_denom=1>

H_0_b3 = '(RD = 0)'
W_test = lm_KM_dyn.wald_test(H_0_b3, use_f = False)
W_test

<class 'statsmodels.stats.contrast.ContrastResults'>
<Wald test (chi2): statistic=[[11.90371726]], p-value=0.0005602203401182512, df_denom=1>

Caso multivariante (especificación VAR): Evolución temporal de tipos de cambio

En esta aplicación se estimará un modelo vectorial autorregresivo para explicar la evolución en el tiempo de las variaciones diarias en los tipos de cambio respecto al dólar de tres de las principales monedas al nivel mundial, el euro, el yen y la libra.

El vector \(\mathbf{z}\) que se usará para el estimar el modelo \(\mathbf{z}_{t} = \mathbf{A}_{0} + \mathbf{A}_{1}\mathbf{z}_{t - 1} + \mathbf{A}_{2}\mathbf{z}_{t - 2} + \ldots\mathbf{A}_{p}\mathbf{z}_{t - p} + \mathbf{u}_{t}\) está compuesto en este caso por los cambios porcentuales diarios de los tipos de cambio, \(\mathbf{z} = (rEurUsd, rUsdJpy, rGbpUsd)'\), donde \(EurUsd\), \(UsdJpy\) y \(GbpUsd\) son, respectivamente, los tipos de cambio de cierre diarios respecto al dólar del euro, el yen y la libra.

Código R

# Lectura de librerías
library(tidyverse)
library(imputeTS)
library(FinTS)
library(vars)
library(rmgarch)
# Lectura de datos
# Descarga de las series de tipos de cambio EURUSD, USDJPY y GBPUSD
# de la página https://finance.yahoo.com/currencies
library(tseries)
EURUSD <- get.hist.quote("EURUSD=X", start = "2003-12-01", end = "2023-12-31")

time series ends   2023-12-29

USDJPY <- get.hist.quote("JPY=X", start = "2003-12-01", end = "2023-12-31")

time series ends   2023-12-29

GBPUSD <- get.hist.quote("GBPUSD=X", start = "2003-12-01", end = "2023-12-31")

time series ends   2023-12-29

XRATES_0 <- merge(EURUSD$Close, USDJPY$Close, GBPUSD$Close)
class(XRATES_0)

[1] "zoo"

autoplot(XRATES_0) + facet_free()

# Análisis de datos perdidos (missing data: NA)
statsNA(EURUSD$Close)

[1] "Length of time series:"
[1] 5240
[1] "-------------------------"
[1] "Number of Missing Values:"
[1] 29
[1] "-------------------------"
[1] "Percentage of Missing Values:"
[1] "0.553%"
[1] "-------------------------"
[1] "Number of Gaps:"
[1] 15
[1] "-------------------------"
[1] "Average Gap Size:"
[1] 1.933333
[1] "-------------------------"
[1] "Stats for Bins"
[1] "  Bin 1 (1310 values from 1 to 1310) :      24 NAs (1.83%)"
[1] "  Bin 2 (1310 values from 1311 to 2620) :      2 NAs (0.153%)"
[1] "  Bin 3 (1310 values from 2621 to 3930) :      2 NAs (0.153%)"
[1] "  Bin 4 (1310 values from 3931 to 5240) :      1 NAs (0.0763%)"
[1] "-------------------------"
[1] "Longest NA gap (series of consecutive NAs)"
[1] "11 in a row"
[1] "-------------------------"
[1] "Most frequent gap size (series of consecutive NA series)"
[1] "1 NA in a row (occurring 13 times)"
[1] "-------------------------"
[1] "Gap size accounting for most NAs"
[1] "1 NA in a row (occurring 13 times, making up for overall 13 NAs)"
[1] "-------------------------"
[1] "Overview NA series"
[1] "  1 NA in a row: 13 times"
[1] "  5 NA in a row: 1 times"
[1] "  11 NA in a row: 1 times"

ggplot_na_distribution(EURUSD$Close)

statsNA(USDJPY$Close)

[1] "Length of time series:"
[1] 5240
[1] "-------------------------"
[1] "Number of Missing Values:"
[1] 29
[1] "-------------------------"
[1] "Percentage of Missing Values:"
[1] "0.553%"
[1] "-------------------------"
[1] "Number of Gaps:"
[1] 15
[1] "-------------------------"
[1] "Average Gap Size:"
[1] 1.933333
[1] "-------------------------"
[1] "Stats for Bins"
[1] "  Bin 1 (1310 values from 1 to 1310) :      24 NAs (1.83%)"
[1] "  Bin 2 (1310 values from 1311 to 2620) :      2 NAs (0.153%)"
[1] "  Bin 3 (1310 values from 2621 to 3930) :      2 NAs (0.153%)"
[1] "  Bin 4 (1310 values from 3931 to 5240) :      1 NAs (0.0763%)"
[1] "-------------------------"
[1] "Longest NA gap (series of consecutive NAs)"
[1] "11 in a row"
[1] "-------------------------"
[1] "Most frequent gap size (series of consecutive NA series)"
[1] "1 NA in a row (occurring 13 times)"
[1] "-------------------------"
[1] "Gap size accounting for most NAs"
[1] "1 NA in a row (occurring 13 times, making up for overall 13 NAs)"
[1] "-------------------------"
[1] "Overview NA series"
[1] "  1 NA in a row: 13 times"
[1] "  5 NA in a row: 1 times"
[1] "  11 NA in a row: 1 times"

ggplot_na_distribution(USDJPY$Close)

statsNA(GBPUSD$Close)

[1] "Length of time series:"
[1] 5240
[1] "-------------------------"
[1] "Number of Missing Values:"
[1] 17
[1] "-------------------------"
[1] "Percentage of Missing Values:"
[1] "0.324%"
[1] "-------------------------"
[1] "Number of Gaps:"
[1] 14
[1] "-------------------------"
[1] "Average Gap Size:"
[1] 1.214286
[1] "-------------------------"
[1] "Stats for Bins"
[1] "  Bin 1 (1310 values from 1 to 1310) :      12 NAs (0.916%)"
[1] "  Bin 2 (1310 values from 1311 to 2620) :      2 NAs (0.153%)"
[1] "  Bin 3 (1310 values from 2621 to 3930) :      2 NAs (0.153%)"
[1] "  Bin 4 (1310 values from 3931 to 5240) :      1 NAs (0.0763%)"
[1] "-------------------------"
[1] "Longest NA gap (series of consecutive NAs)"
[1] "3 in a row"
[1] "-------------------------"
[1] "Most frequent gap size (series of consecutive NA series)"
[1] "1 NA in a row (occurring 12 times)"
[1] "-------------------------"
[1] "Gap size accounting for most NAs"
[1] "1 NA in a row (occurring 12 times, making up for overall 12 NAs)"
[1] "-------------------------"
[1] "Overview NA series"
[1] "  1 NA in a row: 12 times"
[1] "  2 NA in a row: 1 times"
[1] "  3 NA in a row: 1 times"

ggplot_na_distribution(GBPUSD$Close)

# Solución a la presencia de valores ausentes
# Reemplazo de NAs por valores interpolados linealmente
# (la opción na.spline usa splines cúbicos para la interpolación)
EurUsd <- na.approx(EURUSD$Close)
UsdJpy <- na.approx(USDJPY$Close)
GbpUsd <- na.approx(GBPUSD$Close)
XRATES <- merge(EurUsd, UsdJpy, GbpUsd)
autoplot(XRATES) + facet_free()

# NOTA: Se guardan los datos para el análisis en Python
# write.zoo(XRATES, "TIPOS_CAMB.CSV", index.name = "Date")
# Cálculo de las variables del modelo VAR:
# cambios porcentuales diarios (aproximación dif-log)
rEurUsd <- 100 * diff(log(EurUsd))
autoplot(rEurUsd)

rUsdJpy <- 100 * diff(log(UsdJpy))
autoplot(rUsdJpy)

rGbpUsd <- 100 * diff(log(GbpUsd))
autoplot(rGbpUsd)

# Agrupación de los datos en una matriz para el análisis VAR
VAR_data <- as.matrix(cbind(rEurUsd,rUsdJpy,rGbpUsd))
class(VAR_data)

[1] "matrix" "array" 

plot.ts(VAR_data)

# Modelo VAR para las variables rEurUsd, rUsdJpy y rGbpUsd
# Librería vars (https://cran.r-project.org/web/packages/vars/)
# Elección del retardo óptimo (type = c("const", "trend", "both", "none"))
VARselect(VAR_data,lag.max = 5, type = "const")

$selection
AIC(n)  HQ(n)  SC(n) FPE(n) 
     3      2      2      3 

$criteria
                 1           2           3           4           5
AIC(n) -2.80306349 -2.82666306 -2.82771607 -2.82637410 -2.82671302
HQ(n)  -2.79780201 -2.81745546 -2.81456236 -2.80927427 -2.80566707
SC(n)  -2.78801665 -2.80033109 -2.79009897 -2.77747186 -2.76652565
FPE(n)  0.06062406  0.05921011  0.05914779  0.05922722  0.05920715

# Estimación del modelo VAR(2)
VAR2 <-  VAR(VAR_data, p=2)
summary(VAR2)

VAR Estimation Results:
========================= 
Endogenous variables: rEurUsd, rUsdJpy, rGbpUsd 
Deterministic variables: const 
Sample size: 5237 
Log Likelihood: -14867.524 
Roots of the characteristic polynomial:
0.3972 0.3972 0.1218 0.1218 0.09244 0.07614
Call:
VAR(y = VAR_data, p = 2)


Estimation results for equation rEurUsd: 
======================================== 
rEurUsd = rEurUsd.l1 + rUsdJpy.l1 + rGbpUsd.l1 + rEurUsd.l2 + rUsdJpy.l2 + rGbpUsd.l2 + const 

             Estimate Std. Error t value Pr(>|t|)    
rEurUsd.l1 -0.2645415  0.0163622 -16.168  < 2e-16 ***
rUsdJpy.l1 -0.1086175  0.0135507  -8.016 1.34e-15 ***
rGbpUsd.l1  0.1285755  0.0191449   6.716 2.07e-11 ***
rEurUsd.l2 -0.1032622  0.0164544  -6.276 3.76e-10 ***
rUsdJpy.l2 -0.0862361  0.0134700  -6.402 1.67e-10 ***
rGbpUsd.l2  0.0212043  0.0191176   1.109    0.267    
const      -0.0004788  0.0095229  -0.050    0.960    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 0.689 on 5230 degrees of freedom
Multiple R-Squared: 0.066,  Adjusted R-squared: 0.06492 
F-statistic: 61.59 on 6 and 5230 DF,  p-value: < 2.2e-16 


Estimation results for equation rUsdJpy: 
======================================== 
rUsdJpy = rEurUsd.l1 + rUsdJpy.l1 + rGbpUsd.l1 + rEurUsd.l2 + rUsdJpy.l2 + rGbpUsd.l2 + const 

            Estimate Std. Error t value Pr(>|t|)    
rEurUsd.l1 -0.208080   0.016832 -12.362  < 2e-16 ***
rUsdJpy.l1 -0.179985   0.013940 -12.912  < 2e-16 ***
rGbpUsd.l1  0.085397   0.019694   4.336 1.48e-05 ***
rEurUsd.l2 -0.053506   0.016926  -3.161  0.00158 ** 
rUsdJpy.l2 -0.065954   0.013856  -4.760 1.99e-06 ***
rGbpUsd.l2  0.016067   0.019666   0.817  0.41396    
const       0.006500   0.009796   0.664  0.50701    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 0.7088 on 5230 degrees of freedom
Multiple R-Squared: 0.0635, Adjusted R-squared: 0.06242 
F-statistic:  59.1 on 6 and 5230 DF,  p-value: < 2.2e-16 


Estimation results for equation rGbpUsd: 
======================================== 
rGbpUsd = rEurUsd.l1 + rUsdJpy.l1 + rGbpUsd.l1 + rEurUsd.l2 + rUsdJpy.l2 + rGbpUsd.l2 + const 

            Estimate Std. Error t value Pr(>|t|)  
rEurUsd.l1 -0.015425   0.014215  -1.085   0.2779  
rUsdJpy.l1  0.006271   0.011773   0.533   0.5943  
rGbpUsd.l1  0.034416   0.016633   2.069   0.0386 *
rEurUsd.l2 -0.019181   0.014295  -1.342   0.1797  
rUsdJpy.l2 -0.004401   0.011703  -0.376   0.7069  
rGbpUsd.l2  0.008667   0.016609   0.522   0.6018  
const      -0.005646   0.008273  -0.682   0.4950  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 0.5986 on 5230 degrees of freedom
Multiple R-Squared: 0.001177,   Adjusted R-squared: 3.078e-05 
F-statistic: 1.027 on 6 and 5230 DF,  p-value: 0.4055 



Covariance matrix of residuals:
          rEurUsd   rUsdJpy  rGbpUsd
rEurUsd  0.474785 -0.003627  0.22434
rUsdJpy -0.003627  0.502420 -0.05376
rGbpUsd  0.224338 -0.053761  0.35837

Correlation matrix of residuals:
          rEurUsd   rUsdJpy rGbpUsd
rEurUsd  1.000000 -0.007427  0.5439
rUsdJpy -0.007427  1.000000 -0.1267
rGbpUsd  0.543865 -0.126698  1.0000

# Análisis de causalidad de Granger
causality(VAR2, cause = c("rGbpUsd","rUsdJpy"))$Granger

    Granger causality H0: rUsdJpy rGbpUsd do not Granger-cause rEurUsd

data:  VAR object VAR2
F-Test = 38.143, df1 = 4, df2 = 15690, p-value < 2.2e-16

causality(VAR2, cause = c("rEurUsd","rUsdJpy"))$Granger

    Granger causality H0: rEurUsd rUsdJpy do not Granger-cause rGbpUsd

data:  VAR object VAR2
F-Test = 0.78936, df1 = 4, df2 = 15690, p-value = 0.5319

causality(VAR2, cause = c("rEurUsd","rGbpUsd"))$Granger

    Granger causality H0: rEurUsd rGbpUsd do not Granger-cause rUsdJpy

data:  VAR object VAR2
F-Test = 40.419, df1 = 4, df2 = 15690, p-value < 2.2e-16

# Funciones de respuesta al impulso (IRFs)
irf_VAR2 <- irf(VAR2, n.ahead = 5)
plot(irf_VAR2)

# Funciones de descomposición de la varianza (FEVDs)
fevd_VAR2 <- fevd(VAR2, n.ahead = 5)
plot(fevd_VAR2)

# ANEXO: Modelo VAR con estructura GARCH multivariante (MGARCH))
# Librería rmgarch (https://cran.r-project.org/web/packages/rmgarch/)
# Contraste ARCH para los residuos del modelo VAR(2)
ArchTest(VAR2$varresult$rEurUsd$residuals , lags = 1)

    ARCH LM-test; Null hypothesis: no ARCH effects

data:  VAR2$varresult$rEurUsd$residuals
Chi-squared = 434.21, df = 1, p-value < 2.2e-16

ArchTest(VAR2$varresult$rGbpUsd$residuals , lags = 1)

    ARCH LM-test; Null hypothesis: no ARCH effects

data:  VAR2$varresult$rGbpUsd$residuals
Chi-squared = 54.996, df = 1, p-value = 1.208e-13

ArchTest(VAR2$varresult$rUsdJpy$residuals , lags = 1)

    ARCH LM-test; Null hypothesis: no ARCH effects

data:  VAR2$varresult$rUsdJpy$residuals
Chi-squared = 745.32, df = 1, p-value < 2.2e-16

# Modelo VAR con correlación condicional (GARCH-Copula model)
VAR_GARCH_spec  <-  ugarchspec(
variance.model = list (garchOrder =c(1,1), model = "sGARCH") )
Multi_spec  <-  multispec ( replicate (3 , VAR_GARCH_spec) )
Copula_GARCH_spec  <-  cgarchspec (Multi_spec, VAR = TRUE, lag = 2)
VAR_GARCH_fit <-  cgarchfit ( Copula_GARCH_spec , data = VAR_data )
show(VAR_GARCH_fit)

*-------------------------------------------------*
*                  Copula GARCH Fit               *
*-------------------------------------------------*

Distribution        :  mvnorm
No. of Parameters   :  30
[VAR GARCH CC]      : [21+9+0]
No. of Series       :  3
No. of Observations :  5239
Log-Likelihood      :  -12013.46
Av.Log-Likelihood   :  -2.293 

Optimal Parameters
---------------------------------------------------
                  Estimate  Std. Error  t value Pr(>|t|)
[rEurUsd].omega   0.001137    0.000667   1.7043 0.088327
[rEurUsd].alpha1  0.041876    0.003549  11.7979 0.000000
[rEurUsd].beta1   0.957124    0.001271 753.2873 0.000000
[rUsdJpy].omega   0.003261    0.002539   1.2844 0.199017
[rUsdJpy].alpha1  0.065382    0.019181   3.4087 0.000653
[rUsdJpy].beta1   0.931255    0.023867  39.0192 0.000000
[rGbpUsd].omega   0.004389    0.003261   1.3460 0.178315
[rGbpUsd].alpha1  0.057312    0.027181   2.1085 0.034985
[rGbpUsd].beta1   0.931312    0.033560  27.7507 0.000000

Information Criteria
---------------------
                   
Akaike       4.5896
Bayes        4.6009
Shibata      4.5896
Hannan-Quinn 4.5935


Elapsed time : 4.17373 

# Correlaciones condicionales entre los residuos estandarizados
VAR_GARCH_fit@mfit$Rt

           [,1]       [,2]       [,3]
[1,]  1.0000000 -0.2838862  0.6227347
[2,] -0.2838862  1.0000000 -0.2192316
[3,]  0.6227347 -0.2192316  1.0000000

Código Python

# Lectura de librerías
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import statsmodels.tsa.api as smt
# Lectura de datos
XRATES = pd.read_csv("data/TIPOS_CAMB.csv", delimiter = ' ', index_col = 0)
def LogDiff(x):
    x_diff = 100*np.log(x/x.shift(1))
    x_diff = x_diff.dropna()
    return x_diff
VAR_data = pd.DataFrame({'rEurUsd':LogDiff(XRATES['EurUsd']),
                     'rUsdJpy':LogDiff(XRATES['UsdJpy']),
                     'rGbpUsd':LogDiff(XRATES['GbpUsd'])})
# Modelo VAR 
model = smt.VAR(VAR_data)
VARselect = model.select_order(maxlags=5)
print(VARselect.summary())

 VAR Order Selection (* highlights the minimums) 
=================================================
      AIC         BIC         FPE         HQIC   
-------------------------------------------------
0      -2.688      -2.684     0.06805      -2.686
1      -2.803      -2.788     0.06062      -2.798
2      -2.827     -2.800*     0.05921     -2.817*
3     -2.828*      -2.790    0.05915*      -2.815
4      -2.826      -2.777     0.05923      -2.809
5      -2.827      -2.767     0.05921      -2.806
-------------------------------------------------

VAR2 = model.fit(maxlags=2)
print(VAR2.summary())

  Summary of Regression Results   
==================================
Model:                         VAR
Method:                        OLS
Date:           Sun, 09, Feb, 2025
Time:                     14:21:58
--------------------------------------------------------------------
No. of Equations:         3.00000    BIC:                   -2.80141
Nobs:                     5237.00    HQIC:                  -2.81853
Log likelihood:          -14867.5    FPE:                  0.0591468
AIC:                     -2.82773    Det(Omega_mle):       0.0589102
--------------------------------------------------------------------
Results for equation rEurUsd
=============================================================================
                coefficient       std. error           t-stat            prob
-----------------------------------------------------------------------------
const             -0.000479         0.009523           -0.050           0.960
L1.rEurUsd        -0.264542         0.016362          -16.168           0.000
L1.rUsdJpy        -0.108617         0.013551           -8.016           0.000
L1.rGbpUsd         0.128576         0.019145            6.716           0.000
L2.rEurUsd        -0.103262         0.016454           -6.276           0.000
L2.rUsdJpy        -0.086236         0.013470           -6.402           0.000
L2.rGbpUsd         0.021204         0.019118            1.109           0.267
=============================================================================

Results for equation rUsdJpy
=============================================================================
                coefficient       std. error           t-stat            prob
-----------------------------------------------------------------------------
const              0.006500         0.009796            0.664           0.507
L1.rEurUsd        -0.208080         0.016832          -12.362           0.000
L1.rUsdJpy        -0.179985         0.013940          -12.912           0.000
L1.rGbpUsd         0.085397         0.019694            4.336           0.000
L2.rEurUsd        -0.053506         0.016926           -3.161           0.002
L2.rUsdJpy        -0.065954         0.013856           -4.760           0.000
L2.rGbpUsd         0.016067         0.019666            0.817           0.414
=============================================================================

Results for equation rGbpUsd
=============================================================================
                coefficient       std. error           t-stat            prob
-----------------------------------------------------------------------------
const             -0.005646         0.008273           -0.682           0.495
L1.rEurUsd        -0.015425         0.014215           -1.085           0.278
L1.rUsdJpy         0.006271         0.011773            0.533           0.594
L1.rGbpUsd         0.034416         0.016633            2.069           0.039
L2.rEurUsd        -0.019181         0.014295           -1.342           0.180
L2.rUsdJpy        -0.004401         0.011703           -0.376           0.707
L2.rGbpUsd         0.008667         0.016609            0.522           0.602
=============================================================================

Correlation matrix of residuals
            rEurUsd   rUsdJpy   rGbpUsd
rEurUsd    1.000000 -0.007427  0.543865
rUsdJpy   -0.007427  1.000000 -0.126698
rGbpUsd    0.543865 -0.126698  1.000000

# Funciones de respuesta al impulso y de descomposición de la varianza
irf_VAR2 = VAR2.irf(5)
plt.figure(5)
fig = irf_VAR2.plot()
plt.show()

fevd_VAR2 = VAR2.fevd(5)
plt.figure(6)
fig = fevd_VAR2.plot()
plt.show()