Econometría con Matlab: especificación de estructura ARMA y predicción


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®.

 1.      Una vez cargado la matriz de datos (conteniendo una o varias series, ordenadas por columnas).
 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.

Entradas populares de este blog

Valores perdidos (NA) en R: identificación y tratamiento (I)

Ejemplos de tablas en Stata

Métodos y técnicas de Investigación Económica