Esta guía muestra cómo obtener especificaciones ARMA para una serie de datos y
como utilizar dichas especificaciones para obtener predicciones a distintos
horizontes de tiempo. Hacerlo para una sola serie, puede ser relativamente fácil,
generalizarlo para una matriz de datos y distintos horizontes de tiempo es más
interesante. El mismo se puede aplicar
tanto para una serie individual como para una matriz de datos. Los códigos siguiente
se corresponden a modificaciones de códigos kevinsheppard® y de los ejemplos cargados
de Matwork®.
2. Identificar estructura ARMA
usando el criterio AICs
2 2.1. En la primera parte del código se
especifican el numero max de retardo a testear, las opciones del "optimizador" y las matriz (AICs ) donde se guardaran los el Akaike
de cada estimación así como las matrices que guardaran el orden ARMA de cada
serie.
%% Modelizacion de
la media
MAX_AR = 4; MAX_MA = 4;
models = cell(MAX_AR+1,MAX_MA+1);
%
Options to suppress display. ARMA
models are estimated using lsqnonlin
options = optimset('lsqnonlin');
options.MaxIter = 1000;
options.Display='none';
%
Setup arrays to hold values
AICs = zeros(MAX_AR,MAX_MA);
%
Get the effective T since MAX_AR will be held back
T = length(data(:,1)) -
MAX_AR;
orden_ar = [];
orden_ma = [];
|
2.2. Buscando la especificación ARMA
que minimiza el criterio usado
En el
siguiente código se recorre la matriz de variables (llamada dato) con el for
var, el for i va cambiando la especificación AR y j
la MA. Luego
se presenta la especificación testeada (disp.), se estima el modelo ARMA
(armaxfilter), después se calcula el estadístico AICs, estos se guardan en la matriz AICs que se crea en el bloque
anterior. Por último se utiliza find para buscar cual especificación arrojo el menor
AICs y las matrices orden_ar y orden_ama permiten guardar la especificación
ARMA que se corresponde con dicho estadístico.
for var=1:length(data(1,1:end));
n=0;
for i=0:MAX_AR
for j=0:MAX_MA
if i==0 & j==0;
disp('Aqui da error')
else
disp(['AR: ' num2str(i) ' MA: ' num2str(j)])
% The MAX_AR term below enforces
the holdbac, which is requires
% when comparing AR models with
different lag lengths. Failing to
% use this will produce log-likelihods base don
different number of
% observations.
[p, ser,
ll]=armaxfilter(data(:,var),1,i,j,[],[]);
% Store other values for use later
LLs(i+1,j+1) = ll;
AICs(i+1,j+1) = ll - 2 *
(1+i+j);
BICs(i+1,j+1) = ll
- log(T) * (1+i+j);
end
end
end
display('-----------------------------
Variable: '), var
[BIC_AR,BIC_MA] = find(BICs==max(max(BICs)));
disp('Orden del modelo ARMA: Criterio AICs')
disp([BIC_AR
- 1 BIC_MA - 1])
orden_ar=[orden_ar, BIC_AR-1];
orden_ma=[
orden_ma, BIC_MA-1];
end
|
3. Usar la estructura ARMA para
predecir
Ahora
se puede utilizar esta especiación para obtener predicciones de las series.
-
Nomvar: Crea nombres de las series, para
identificar los gráficos (habrán tantos nombres como columna tenga nuestra
matriz data)
-
For
variable, permite recorrer las series ordenadas en la matriz data(length(data(1,:)) indica que debe ir hasta la ultima columna)
-
arima
Permite estimar el modelo ARIMA llamado Mdl donde se utilizan las
especificaciones ARMA de cada serie, que encontramos en los bloques de códigos anteriores.
(Una aclaración importante, es que se supones series integradas por donde el
orden de esta es cero para todos los casos, esto dependerá de los datos con que
trabajes)
-
ventana
indicas el horizonte de tiempo que quieres que tu modelo prediga (4
indica que el primer grafico mostrara las predicciones del ARIMA para cuatro)
-
no_proyec, Va recorriendo las ventanas de
predicciones
-
estimate , Estima el ARIMA en el punto
del tiempo determinado, es como si estimáramos el modelo el dia en que inicia
la predicción sin conocer lo que sucederá.
-
forecast , Permite obtener las
predicciones del modelo
-
Por
ultimo se construye el gráfico
%% Forescat mean response
nomvar = {'PIB', 'Net Change-Off', 'Desempleo', 'Inflación'...
'Tipo a
Largo', 'Tipo a
Corto'};
for variable=1:length(data(1,:));
Y=data(:,variable);
[nobs, nvar]=size(Y);
Mdl =
arima(orden_ar(variable),0,orden_ma(variable));
ventana= [4, 8, 31];
for no_proyec = 1:length(ventana);
EstMdl =
estimate(Mdl,Y(1:end-ventana(no_proyec),:));
[YF YMSE] =
forecast(EstMdl,ventana(no_proyec),'Y0',Y(1:end-ventana(no_proyec)));
figure(variable)
subplot(1,length(ventana),no_proyec)
h1 = plot(Y,'Color',[.7,.7,.7]);
hold on
inicio=nobs-(ventana(no_proyec)-1);
h2 = plot(inicio:nobs,YF,'b','LineWidth',2);
h3 = plot(inicio:nobs,YF + 1.96*sqrt(YMSE),'r:',...
'LineWidth',2);
plot(inicio:nobs,YF - 1.96*sqrt(YMSE),'r:','LineWidth',2);
legend([h1 h2 h3],'Observed','Forecast',...
'95% Confidence Interval','Location','NorthWest');
legend('boxoff')
title({[num2str(ventana(no_proyec)) '-Period Forecasts and Approximate
95%, for: ',num2str(nomvar{variable})]})
hold off
end
end
|
Tomando
datos trimestrales del Net Change-off de USA y una especificación ARIMA(1,0,4),
se obtiene el siguiente gráfico, en el mismo se muestra la predicción del modelo en tres puntos del tiempo, el ultimo año, los ultimo dos años, y lo que predecirla antes de la crisis.:
Predicciones ARIMA(1,0,4) de probabilidades de
Default
(Datos trimestrales
1984-2014)
Fuente: Elaboración propia con datos de FDIC
Una aclaración al final, es que la especificación ARMA de los primeros bloques es del conjunto completo de datos, pudiendo ser esta cambiante en el tiempo. Habría que verificar si en 2006, por ejemplo la especificación ARMA hubiese sido la que se obtiene usando el periodo completo.
