Forecast: problemas con omitidos (no estimados) en el output

[atorre]

Hola MMS,
en el objeto Forecast adjunto, podéis encontrar un ejemplo de la previsión de un modelo con un único submodelo que tiene un valor omitido, no estimado, en el output. Esto provoca que la previsión no se ejecute.

Esto puede ocurrir cuando se ejecuta la previsión de un modelo después de actualizar los datos y no se ha estimado (normalmente, por falta de tiempo). Es necesario adaptar el mecanismo de previsión para tener en cuenta esto y solucionarlo para poder ejecutar la previsión sin tener en cuenta esos datos desconocidos.

Una posible solución puede verse el código que sigue a continuación rescatado de SADD. La función Forecast ejecuta una previsión y utiliza las funciones ActRes y OmitFilterTOL.

//////////////////////////////////////////////////////////////////////////////
   Serie InterRecur(Set impulsos, Serie resi, Polyn xpiw)
// PURPOSE: Reduce a cero los residuos correspondientes a las fechas de 
//          <impulsos> trasmitiendo el valor del residuo segun la funcion de 
//          respuesta <xpiw>
//////////////////////////////////////////////////////////////////////////////
{
  Date iesima = impulsos[1];
  Real valor  = SerDat(resi, iesima);
  Serie pulEx = xpiw:(valor * Pulse(iesima, Dating(resi)));
  Serie res1  = resi - pulEx; 
  Set imp1    = impulsos - SetOfDate(iesima);
  If(GT(Card(imp1),0), InterRecur(imp1,res1,xpiw), res1)
};

//////////////////////////////////////////////////////////////////////////////
Serie OmitFilterTOL(Ration piW, Serie res, Serie ind)
//////////////////////////////////////////////////////////////////////////////
{
  Polyn xpiW     = Expand(piW,CountS(ind));
  Set pulsosF    = Dates(SerTms(ind),First(ind),Last(ind));
  Real cardpul   = Card(pulsosF);
  If(cardpul, InterRecur(pulsosF,res,xpiW), res)
};
//////////////////////////////////////////////////////////////////////////////
PutDescription(
"funcion que devuelve la serie de correccion de omitidos en la serie de
residuos actualizados.",
OmitFilterTOL);
//////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////
Serie ActRes(Serie serRuido, Serie serRes, Real p, Real q, Ration piw)
//////////////////////////////////////////////////////////////////////////////
{
  Date  lastRui  = Last(serRuido);
  Date  lastRes  = Last(serRes);
  Real  lag      = DateDif(Dating(serRuido), lastRes, lastRui);
  Real  m        = p + lag;
  Serie aIni     = If(q,LastM(serRes,q),CalInd(W,Dating(serRuido)));
  Serie noiseIni = LastM(serRuido,m);
  If(q,DifEq(piw, noiseIni, aIni), Numerator(piw):noiseIni)
};
//////////////////////////////////////////////////////////////////////////////
PutDescription(
"funcion que actualiza los residuos, a partir de los ordenAR ultimos valores 
del ruido y los ordenMA ultimos valores de los residuos. Sus argumentos son: 
  - serRuido, la serie de ruido del modelo estimado
  - serRes, la serie de residuos que se desea actualizar.",
ActRes);
//////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////
Set Forecast (Serie original,   // Serie original
              Set   boxCox,     // Transformacion Box-Cox
              Set   arima,      // Arima Estimado
              Serie noise,      // Noise Estimado
              Serie filter,     // Filtro Extendido
              Serie filterVar,  // SetSum(ForErr(Inputs Estoc)^2)
              Serie residuals,  // Residuos Estimados
              Real  numFor,     // Numero fde previsiones
              Date  origen,     // Origen de prevision (ultimo dato conocido)
              Real  alfa,       // 2*(1-DistNormal(1))
              Serie cierres)    // Serie de cierres de la serie si los tiene
//////////////////////////////////////////////////////////////////////////////
{
  Text  dtnTxt   = "TimeSet "+ DatingName(original) +"; ";

  Polyn   ARI    = ARIMAGetARI(arima);
                                 //arima[2]*arima[4];
  Polyn   MA     = ARIMAGetMA (arima);//arima[3];//
  Real    degARI = Degree(ARI);
  Real    degMA  = Degree(MA);
  Ration  piw    = ARI/MA;
  Ration  psiw   = MA/ARI;

  Date    fsFor  = Succ(origen, Eval(dtnTxt), 1);      // primera prevision 
  Date    lsFor  = Succ(origen, Eval(dtnTxt), numFor); // ultima  prevision
  Real    isEst  = Last(noise) == origen;  // El origen de prevision coincide 
                                          // con el ultimo dato de estimacion
  ////////////////////////////////////////////////////////////////////////////
  // Output  y Noise extendido y Residuos actualizados
  ////////////////////////////////////////////////////////////////////////////
  Serie outputOmit = BoxCoxTransf(original, boxCox);
  Serie indOmi     =  // omitidos no estimados
    If(isEst, CalInd(W,Eval(dtnTxt)),
              SubSer(IsUnknown(outputOmit), 
                      Succ(Last(noise),Eval(dtnTxt),1), 
                      origen));

  Serie output   = IfSer(IsUnknown(outputOmit), 0, outputOmit);
  Serie filter2For = SubSer(filter, fsFor, lsFor);
  Serie noiseExt = noise >> (output - filter);

  Date  lastAvailable    = Min(Last(noiseExt), origen); // ultimo noise
  Date  lastAvailableRes = Min(Last(residuals),origen); // ultimo residuo

  Serie noiseNoFil = SubSer(noiseExt, First(noiseExt), lastAvailable); 
  Serie resEst     = SubSer(residuals, First(residuals), lastAvailableRes); 
  Real  doAct      = Or((Last(noise) < lastAvailable), 
                         Last(residuals) < lastAvailableRes);// Ojo con esto
  Serie actRes     = If(Not(doAct), resEst, // actualizacion de residuos
                            ActRes(noiseNoFil, resEst, degARI, degMA, piw));

  ////////////////////////////////////////////////////////////////////////////
  // Residuos y Ruido Filtrado
  ////////////////////////////////////////////////////////////////////////////
  Date  lsDat     = Last(noiseNoFil);
  Serie filtRes   = If(Not(doAct), actRes, OmitFilterTOL(piw, actRes, indOmi));
  Serie resActF   = actRes - filtRes;
  Serie resActIni = SubSer(CalInd(W,Eval(dtnTxt)),
                          Succ(First(resActF), Eval(dtnTxt), -degMA),
                          Succ(First(resActF), Eval(dtnTxt), -1));
  Serie filNoi   = DifEq(psiw, resActIni<<resActF);
  Serie cierre   = SubSer(CalInd(W,Eval(dtnTxt)), 
                   Succ(lsDat, Eval(dtnTxt), -degARI), lsDat) 
                   << filNoi;

  Serie hisRes = If(degMA,
                    SubSer(filtRes, Succ(lsDat,Eval(dtnTxt),-degMA), lsDat),
                    CalInd(W,Dating(filtRes)));

  Serie hisNoi = SubSer((noiseNoFil - cierre), 
                  Succ(lsDat,Eval(dtnTxt),-degARI), lsDat);

  ////////////////////////////////////////////////////////////////////////////
  // Prevision del noise
  ////////////////////////////////////////////////////////////////////////////
  Serie zero      = SubSer(CalInd(W, Eval(dtnTxt)), fsFor, lsFor);
  Serie allRes    = If(degMA,(hisRes>> zero),zero);
  Serie prevNoise = SubSer(DifEq(psiw, allRes, hisNoi), fsFor, lsFor);
  
  Serie prevTrans = filter2For + prevNoise;

  ////////////////////////////////////////////////////////////////////////////
  // Error de Prevision
  ////////////////////////////////////////////////////////////////////////////
  Real  sigma     = StDsS(hisRes);
  Serie forErrOut = ErrorForecast(numFor, sigma, origen, psiw, Eval(dtnTxt));
 // Set forErrorSet = SetOfSet(ErrorForecast2(numFor, sigma, origen, psiw, Eval(dtnTxt)));
  Serie forErr    = SubSer(SqRt(forErrOut^2 + filterVar), fsFor, lsFor);

  ////////////////////////////////////////////////////////////////////////////
  // Bandas
  ////////////////////////////////////////////////////////////////////////////
  Real  numSig  = DistNormalInv(1-(alfa/2));
  Serie tSup    = prevTrans + (numSig*forErr);
  Serie tInf    = prevTrans - (numSig*forErr);

  ////////////////////////////////////////////////////////////////////////////
  // Salida de la funcion
  ////////////////////////////////////////////////////////////////////////////
  Serie forecast      = BoxCoxInvTransf(prevTrans, boxCox)*Not(cierres);
  Serie viewForecast  = original << forecast;  
//  Serie viewForecast  = 
//   SubSer(original<<forecast, Last(original), Last(forecast));  
  Serie upperBand     = BoxCoxInvTransf(tSup, boxCox);
  Serie lowerBand     = BoxCoxInvTransf(tInf, boxCox);
  Set  previsionResults = 
  [[ original, viewForecast, upperBand, forecast, lowerBand,
     prevNoise, filter2For, prevTrans, forErr]];

  Set realSet   = SetOfReal(numFor);
  Set aRIMASet  = SetOfPolyn(ARI, MA);
  Set dateSet   = SetOfDate(origen, lsFor, lastAvailable, lastAvailableRes);
  Set serieSet  = SetOfSerie
  (
    outputOmit, indOmi, output, noiseExt, noiseNoFil, resEst, 
    actRes, zero, filtRes, resActF, resActIni, filNoi, cierre,
    allRes, hisNoi, forErrOut, forErr, tSup,tInf
  );
  Set diagnosisSet = SetOfSet(aRIMASet, dateSet, serieSet, realSet);// forErrorSet
  previsionResults << [[ diagnosisSet ]]
};