<jupyter_start><jupyter_text>Modelo VAR<jupyter_code>import pandas as pd
import numpy as np
import random
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn; seaborn.set()
plt.rcParams['figure.figsize'] = [10, 8]  #por defecto las figuras se displayan 8 pulgadas de ancho y 6 de alto<jupyter_output><empty_output><jupyter_text>VARAbajo se presenta un ejemplo de construcción de modelo VAR para el dataset macrodata provisto en statsmodels.api.   Dicho dataset tiene 14 variables económicas de USA a una resolución trimestral recolectadod por la FRED. Para este ejemplo se trabajará solo con 3 variables:    * realgdp:   Producto interno bruto real de EEUU (Bil. USD encadenados de 2005, tasa anual desestacionalizada)    * realcons:  Gastos reales de consumo personal  (Bil. USD encadenados de 2005, tasa anual desestacionalizada)    * realinv:   Inversión interna privada bruta real (Bil. USD encadenados de 2005, tasa anual desestacionalizada)    Mayor informacion puede ser encontrada en https://www.statsmodels.org/dev/datasets/generated/macrodata.html<jupyter_code># Librerias a usar
import numpy as np
import pandas
import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.tsa.api import VAR
from statsmodels.tsa.base.datetools import dates_from_str

# Lee la data
mdata = sm.datasets.macrodata.load_pandas().data

# prepara el indice temporal
dates = mdata[['year', 'quarter']].astype(int).astype(str)
quarterly = dates["year"] + "Q" + dates["quarter"]
quarterly = dates_from_str(quarterly)

# extrae las 3 variables que se trabajará y coloca el indice temporal
mdata = mdata[['realgdp','realcons','realinv']]
mdata.index = pandas.DatetimeIndex(quarterly)

# Retira la tendencia haciendo diferencias y aplica transformacion logaritmica para estacionarizar la varianza 
data = np.log(mdata).diff().dropna()

# Instancia el modelo VAR
model = VAR(data)
# plotea las variables originales (sin transformar)
mdata.plot()
plt.show()
# Plotea la data transformada (sin tendencias y estacionarizada en varianza)
data.plot()
plt.show()
# Computa varios VAR  (hasta 10 lags) e informa las configuraciones con mejores valores de criterios de Informacion  
print(model.select_order(10))
# Ajusta un modelo VAR con p=1 (fue el lag con mejores criterios AIC, BIC)
model_fit = model.fit(1)
# Sumariza el modelo ajustado
model_fit.summary()

# Plotea autocorrelación y crosscorrelacion de los residuos (para ver si vale la pena un VARMA) 
model_fit.plot_acorr(resid=False)
# predice el siguiente punto temporal en las 3 variables (hay que proporcionarle la data histórica y el numero de steps)
y_predicted = model_fit.forecast(data.values[-1:], 1)
data.values[-2:]
y_predicted
# predice los siguientes 50 puntos (es predicción encadenada, no es walk-forwad ) 
model_fit.plot_forecast(10)<jupyter_output><empty_output><jupyter_text>Se busca- Implementar la evaluacion en walk-forward del modelo VAR en el dataset macrodata  usando los 50 ultimos puntos como test. Calcular el rmse de cada variable- Comparar con modelos AR ajustados individualmente en cada una de las variables individualmente      VAR walk-forward<jupyter_code>import warnings
from sklearn.metrics import mean_squared_error   # para calcular el error cuadratico medio
from math import sqrt
warnings.filterwarnings("ignore")   # para evitar que muestre constantemente advertencias

#50 registros para el test
train,test = data[0:-50],data[-50:]

historia=train
predicciones = pd.DataFrame()
# evaluación walk-forward
for i in range(len(test)):
# Instancia el modelo VAR con la historia
  modelo_wf = VAR(historia)
# Computa varios VAR  (hasta 10 lags) e informa las configuraciones con mejores valores de criterios de Informacion  
#  print(modelo_wf.select_order(10))
# Ajusta un modelo VAR con p=1 (fue el lag con mejores criterios AIC, BIC)
  modelo_wf_fit = modelo_wf.fit(1)
# predice el siguiente punto temporal en las 3 variables (hay que proporcionarle la data histórica (solo un lag, en este caso) y el numero de steps (=1))
  y_predicted = model_fit.forecast(historia.values[-1:], 1)  
  predicciones = predicciones.append(pd.DataFrame(y_predicted), ignore_index = True)
  historia=historia.append(test.iloc[i])

# evalua el RMSE para las 50 predicciones por cada variable
for j in range(0,3):
 rmse = sqrt(mean_squared_error(test.iloc[:,j], predicciones[j]))
 print('RMSE walk-forward en test: %.3f' % rmse,'en variable',test.columns[j])<jupyter_output>RMSE walk-forward en test: 0.005 en variable realgdp
RMSE walk-forward en test: 0.005 en variable realcons
RMSE walk-forward en test: 0.031 en variable realinv<jupyter_text>Modelo AR - Variable realgdp<jupyter_code>from statsmodels.tsa.ar_model import AR
# entrena un modelo AR
ar_model = AR(train.iloc[:,0])
ar_model_fit =  ar_model.fit()

# imprime el orden encontrado (p) del modelo y sus coeficientes
p = ar_model_fit.k_ar
coef = ar_model_fit.params
print('Orden p del modelo AR: %s' % p)
print('Coeficientes del modelo AR: %s' % ar_model_fit.params)
# Funcion para predecir el siguiente valor con los coeficientes de un modelo AR y una lista de observaciones históricas
def predict(coef, history):
    y_predicted = coef[0]
    for i in range(1, len(coef)):
        y_predicted += coef[i] * history[-i]
    return y_predicted

# hace las predicciones en test
history = pd.Series(train.iloc[:,0]).tolist()

predictions = list()
for t in range(len(test)):
    y_predicted = predict(coef, history)
    predictions.append(y_predicted)
    history.append( pd.Series(test.iloc[t,0]) )  # se añade la ultima observacion real al histórico (al estilo walk-forward)     
        
rmse = sqrt(mean_squared_error(test.iloc[:,0], predictions))
print('RMSE en test: %.3f' % rmse)<jupyter_output>RMSE en test: 0.006<jupyter_text>Modelo AR - Variable realinv<jupyter_code>from statsmodels.tsa.ar_model import AR
# entrena un modelo AR
ar_model = AR(train.iloc[:,2])
ar_model_fit =  ar_model.fit()

# imprime el orden encontrado (p) del modelo y sus coeficientes
p = ar_model_fit.k_ar
coef = ar_model_fit.params
print('Orden p del modelo AR: %s' % p)
print('Coeficientes del modelo AR: %s' % ar_model_fit.params)
# Funcion para predecir el siguiente valor con los coeficientes de un modelo AR y una lista de observaciones históricas
def predict(coef, history):
    y_predicted = coef[0]
    for i in range(1, len(coef)):
        y_predicted += coef[i] * history[-i]
    return y_predicted

# hace las predicciones en test
history = pd.Series(train.iloc[:,1]).tolist()

predictions = list()
for t in range(len(test)):
    y_predicted = predict(coef, history)
    predictions.append(y_predicted)
    history.append( pd.Series(test.iloc[t,1]) )  # se añade la ultima observacion real al histórico (al estilo walk-forward) 
        
rmse = sqrt(mean_squared_error(test.iloc[:,1], predictions))
print('RMSE en test: %.3f' % rmse)<jupyter_output>RMSE en test: 0.005<jupyter_text>Modelo AR - Variable realinv<jupyter_code>from statsmodels.tsa.ar_model import AR
# entrena un modelo AR
ar_model = AR(train.iloc[:,2])
ar_model_fit =  ar_model.fit()

# imprime el orden encontrado (p) del modelo y sus coeficientes
p = ar_model_fit.k_ar
coef = ar_model_fit.params
print('Orden p del modelo AR: %s' % p)
print('Coeficientes del modelo AR: %s' % ar_model_fit.params)
# Funcion para predecir el siguiente valor con los coeficientes de un modelo AR y una lista de observaciones históricas
def predict(coef, history):
    y_predicted = coef[0]
    for i in range(1, len(coef)):
        y_predicted += coef[i] * history[-i]
    return y_predicted

# hace las predicciones en test
history = pd.Series(train.iloc[:,2]).tolist()

predictions = list()
for t in range(len(test)):
    y_predicted = predict(coef, history)
    predictions.append(y_predicted)
    history.append( pd.Series(test.iloc[t,2]) )  # se añade la ultima observacion real al histórico (al estilo walk-forward) 
        
rmse = sqrt(mean_squared_error(test.iloc[:,2], predictions))
print('RMSE en test: %.3f' % rmse)<jupyter_output>RMSE en test: 0.036
{'fraction_non_alphanumeric': 0.09347420758234928, 'fraction_numerical': 0.012927284027346178, 'mean_word_length': 5.031484257871065, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 21, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '148625', 'license': 'MIT', 'n_tokens_mistral': 2873, 'n_tokens_neox': 2800, 'n_words': 872, 'path': 'time-series/ModeloVAR.ipynb', 'repo_name': 'rasecotineb/datascience', 'size': '8045'}