Ticket #904: _decision.tol

//////////////////////////////////////////////////////////////////////////////
// FILE    : _decision.tol
// AUTHOR : Luis Miguel Pérez Liberal (lmperez@bayesinf.com)
// PURPOSE : Standard functions for diagnosis procedures
//////////////////////////////////////////////////////////////////////////////


//////////////////////////////////////////////////////////////////////////////
Serie ErrorForecast.Clasico(Real lengthFor, Serie residuals, Ratio psiw)
//////////////////////////////////////////////////////////////////////////////
{
  TimeSet dtn = Dating(residuals);
  Date origen = Last(residuals);
  Real sigma  = StDsS(residuals);

  Real  sigma2  = sigma^2;
  Date  fsFor   = Succ(origen, dtn, 1);
  Date  lsFor   = Succ(origen, dtn, lengthFor+1);
  Polyn psi     = Expand(psiw, lengthFor);
  Serie zero    = CalInd(W,dtn);
  Serie psiSer  = SubSer(psi:(zero + Pulse(origen,dtn)), origen, lsFor);
  Serie psiSer2 = psiSer^2;
  Serie error2  = sigma2 * DifEq(1/(1-B), psiSer2, 0);

  Serie error   = SqRt(error2)
};
//////////////////////////////////////////////////////////////////////////////
PutDescription("Devuelve la sigma en prevision de una serie de residuos.
Ejemplo:

Serie ErrorForecast.Clasico(
  90,
  GenResiduals.Heteroscedastic(500),
  Ratio (1-0.75*B^7)/((1-B^7)*(1-0.10*B))
);
",
ErrorForecast.Clasico);
//////////////////////////////////////////////////////////////////////////////


//////////////////////////////////////////////////////////////////////////////
Serie ErrorForecast.Het.Luis(Real lengthFor, Serie residuals, Ratio psiw)
//////////////////////////////////////////////////////////////////////////////
{
  TimeSet dtn = Dating(residuals);

  Set sTms = {For(1,7,Text (Real k){ 
    Text If(CountS(DatCh(residuals,WD(k),FirstS)),"WD("+IntText(k)+")","")
  }) - [[ "" ]]};

  Date origen   = Last(residuals);
  Date  fsFor   = Succ(origen, dtn, 1);
  Date  lsFor   = Succ(origen, dtn, lengthFor+1);
  Polyn psi     = Expand(psiw, lengthFor);
  Serie zero    = CalInd(W,dtn);
  Serie psiSer  = SubSer(psi:(zero + Pulse(origen,dtn)), origen, lsFor);
  Serie psiSer2 = psiSer^2;
  Serie difEQ   = DifEq(1/(1-B), psiSer2, 0);

  Serie error2  = SetSumC(EvalSet(sTms,Serie (Text tms)
  {
    TimeSet tms_WD   = Eval(tms);
    Real  sigma_WD   = StDsS(DatCh(residuals,tms_WD,FirstS));
    Real  sigma2_WD  = sigma_WD^2;
    Serie psiSer2_WD = difEQ*sigma2_WD*CalInd(tms_WD,dtn);
    Serie psiSer2_WD
  }));

  Serie DatCh(SqRt(SubSer(error2,fsFor,lsFor)),Dating(residuals),FirstS)
};
//////////////////////////////////////////////////////////////////////////////
PutDescription("Devuelve la sigma en prevision de una serie de residuos en
fechado Diario considerando la heterocedasticidad de los mismos.
Esta función falla para modelos con parámetros arima regulares grandes.
Ejemplo:

Serie ErrorForecast.Het.Luis(
  90,
  GenResiduals.Heteroscedastic(500),
  Ratio (1-0.75*B^7)/((1-B^7)*(1-0.10*B))
);",
ErrorForecast.Het.Luis);
//////////////////////////////////////////////////////////////////////////////


//////////////////////////////////////////////////////////////////////////////
Serie ErrorForecast.Het.Pepe(Real lengthFor, Serie residuals, Serie noise,
  Set arima)
//////////////////////////////////////////////////////////////////////////////
{
  TimeSet dtn = Dating(residuals);
  Date origen = Last(residuals);
  Date  fsFor = Succ(origen, dtn, 1);
  Date  lsFor = Succ(origen, dtn, lengthFor+1);

  Set   sAr   = arima[2];
  Set   sMa   = arima[3];
  Polyn dif   = SetProd(arima[4]);
  Ratio psiw  = ma/(dif*ar);
  Real maDeg  = Degree(SetProd(sMa));
  Real ariDeg = Degree(SetProd(sAr)*dif);

  Set sTms = {For(1,7,Text (Real k){ 
    Text If(CountS(DatCh(residuals,WD(k),FirstS)),"WD("+IntText(k)+")","")
  }) - [[ "" ]]};

  Serie sigmaK = SetSumC(EvalSet(sTms,Serie (Text tms){
    Serie CalInd(Eval(tms),dtn)*StDsS(DatCh(residuals,Eval(tms),FirstS))
  }));

  Real noDegree = If(GT(maDeg+ariDeg,CountS(noise)),1,0);

  Serie Case(
    EQ(MaxS(sigmaK),0), ErrorForecast.Clasico(lengthFor,residuals,psiw),
    EQ(noDegree,1),     ErrorForecast.Het.Luis(lengthFor,residuals,psiw),
    True,
    {
      Real Show(False,"ALL",True);
      Set estimate = Estimate(
        ModelDef(
          noise/(sigmaK+Not(sigmaK)),
          1,
          0,
          SetMax(arima[1]),
          0,
          dif,
          sAr,
          sMa,
          Empty,
          Empty
        )
      );
      Real Show(True,"ALL");
    
      Polyn p.ar   = SetProd(estimate::Definition::AR);
      Polyn p.ma   = SetProd(estimate::Definition::MA);
      Polyn p.dif  = estimate::Definition::Dif;
      Ratio r.psiw = p.ma/(p.dif*p.ar);
      Serie error  = estimate::Series::Residuals;
  
      Serie error.Prevision = ErrorForecast.Clasico(
        lengthFor,
        error,
        r.psiw
      )/(StDsS(error)+Not(StDsS(error)));
    
      Serie DatCh(
        SubSer(error.Prevision*sigmaK,fsFor,lsFor),
        Dating(residuals),
        FirstS
      )
    }
  )
};
//////////////////////////////////////////////////////////////////////////////
PutDescription("Devuelve la sigma en prevision de una serie de residuos en
fechado Diario considerando la heterocedasticidad de los mismos.
Esta función es capaz de controlar la heterocedasticidad de cualquier modelo
arima; sin embargo, puede darnos problemas si tenemos pocos grados de libertad
y es mucho más costosa computacionalmente.
Ejemplo:

Set sArima  = GetArimaFromLabel(\"P1_7DIF0_1AR1_0MA0_7\");
Polyn p.AR  = SetProd(Traspose(sArima)[2]);
Polyn p.MA  = SetProd(Traspose(sArima)[3]);
Polyn p.DIF = SetProd(Traspose(sArima)[4]);
Polyn ar    = GetRandomPolyn(p.AR);
Polyn ma    = GetRandomPolyn(p.MA);
Ratio psi   = ma/(p.DIF*ar);

Serie residuals = GenResiduals.Heteroscedastic(500);
Serie noise     = DifEq(psi,residuals,100);
  
Serie ErrorForecast.Het.Pepe(
  90,
  residuals,
  noise,
  Traspose(sArima)
);",
ErrorForecast.Het.Pepe);
//////////////////////////////////////////////////////////////////////////////


//////////////////////////////////////////////////////////////////////////////
// END
//////////////////////////////////////////////////////////////////////////////