Ticket #904: _posEstimation.tol

//////////////////////////////////////////////////////////////////////////////
// FILE    : _posEstimation.tol
// AUTHOR : Luis Miguel Pérez Liberal (lmperez@bayesinf.com)
// PURPOSE : Post-Estimation diagnosis functions
//////////////////////////////////////////////////////////////////////////////


//////////////////////////////////////////////////////////////////////////////
Set BPL(Serie serie, Text arima, Real autocorr, Real alpha)
//////////////////////////////////////////////////////////////////////////////
{
  Text label = If(arima=="ARIMA","P1_7DIF0_1AR1_0MA0_7",arima);
  Set  arima = GetArimaFromLabel(label);

  Real rAR = Card(Monomes(ARIMAGetAR(arima))-[[ Polyn 1 ]]);
  Real rMA = Card(Monomes(ARIMAGetMA(arima))-[[ Polyn 1 ]]);
  Real deg = Degree(ARIMAGetDIF(arima));

  Real bpl   = BoxPierceLjung(serie,autocorr,rAR,rMA);
  Real bpl_D = BoxPierceLjung(serie,autocorr*deg,rAR,rMA);

  Real gl  = autocorr     - rAR - rMA;
  Real glD = autocorr*deg - rAR - rMA;

  Real chi  = DistChiInv(1-alpha/2,gl);
  Real chiD = DistChiInv(1-alpha/2,glD);

  Real t_blp   = GT( Abs(bpl), chi );
  Real t_blp_D = GT( Abs(bpl_D),chiD );

  Set Eval("Set _.BPL."+Name(serie)+" = [[ t_blp, t_blp_D ]];")
};
//////////////////////////////////////////////////////////////////////////////
PutDescription("Realiza un contraste de hipotesis BoxPierceLjung sobre las
autocorrelaciones dadas (autocorr) dado un arima (como etiqueta).
Devuelve verdadero si se rechaza la hipotesis y falso si se acepta.
t_blp es el resultado para las autocorrelaciones dadas, y t_blp_D es el
resultado del test para las autocorrelaciones dadas multiplicada por el grado
de la diferencia del arima.",
BPL);
//////////////////////////////////////////////////////////////////////////////


//////////////////////////////////////////////////////////////////////////////
Set Analisis.Residuals(Serie serie,Real nParametros, Real nParametrosArima,
  Real alpha)
//////////////////////////////////////////////////////////////////////////////
{
  Real med = AvrS(serie);
  Real sig = StDsS(serie);
  Real asi = AsimetryS(serie);
  Real kur = KurtosisS(serie);

  Real asi_teo = AsimetryS(Gaussian(med,sig));
  Real kur_teo = KurtosisS(Gaussian(med,sig));

  Real t     = med/sig;
  Real gl    = CountS(serie) - nParametrosArima - nParametros;
  //Real alpha = 0.10;
  Real p_med = DistTInv(1-alpha/2,gl);

  Real test_med = GT( Abs(t),p_med );

  Serie uno  = Not(serie)+Not(Not(serie));
  Serie norm = (serie - med*uno)/StDsS(serie);
  Real  sig1 = SumS(IfSer(GT(Abs(norm),uno*1),1,0))/SumS(uno);
  Real  sig2 = SumS(IfSer(GT(Abs(norm),uno*2),1,0))/SumS(uno);
  Real  sig3 = SumS(IfSer(GT(Abs(norm),uno*3),1,0))/SumS(uno);
  Set   end  = [[ med,sig,asi,asi_teo,kur,kur_teo,test_med,sig1,sig2,sig3 ]];

  Set Eval("Set "+Name(serie)+" = end;")
};
//////////////////////////////////////////////////////////////////////////////
PutDescription("Dada una serie de residuos, devuelve:
Media
Sigma
Asimetria y Asimetría teórica
Kurtosis y Kurtosis teórica
Test media distinta de cero
Numero de residuos fuera de una, dos y tres sigmas.

Precisa que le digamos el alpha como el área fuera de la curva para el test 
de la media",
Analisis.Residuals);
//////////////////////////////////////////////////////////////////////////////


//////////////////////////////////////////////////////////////////////////////
Set Analisis.Rachas(Serie residuos, Real max)
//////////////////////////////////////////////////////////////////////////////
{
  // Crea un polinomio del tipo 1+B+B^2+...+B^k
  Polyn PolAcu_K(Polyn operador,Real k){
    SetSum( For(0,k,Polyn (Real j){ operador^j }) ) 
  };

  // Cuenta las fechas de una serie de rachas, cuenta hacia delante y hacia
  // atras para evitar que se elimine la primera fecha de la racha
  Set DatesRacha(Serie racha)
  {
    Set For(1,max,Set (Real j)
    {
      Serie I   = Not(racha)+Not(Not(racha));
      Serie s1  = IfSer(EQ(PolAcu_K(B,j):racha,I*(j+1)),1,0);
      Serie s2  = IfSer(EQ(PolAcu_K(F,j):racha,I*(j+1)),1,0);
      Set dates = DatesInd(s1+s2)
    })
  };

  // Calcula el numero de ocurrencias de cada racha
  Set DistRacha(Set datesRacha, Serie residuos, Text signo)
  {
    Real med = AvrS(residuos);
    // cojo la sigma de los residuos porque al caso lo que me interesa es la
    // media, y así no tengo que tener en cuenta la transformación
    Real sig = StDsS(residuos);
    Real prb = Case(
      signo=="+", 1-DistNormal(0,med,sig),
      signo=="-", DistNormal(0,med,sig),
      True,"ERROR"
    );

    Text label = Case(
      signo=="+", "pos",
      signo=="-", "neg",
      True,"ERROR"
    );

    Set For(1,Card(datesRacha)-1,Set (Real j)
    {
      Set  dates     = datesRacha[j] - datesRacha[j+1];
      Real observado = Card(dates);
      // No es el número de observaciones, hay que restarle la racha menos 1
      Real longitud  = CountS(residuos)-(j+1-1);
      Real teorico   = Round((prb^(j+1))*longitud,2);
      Set Eval("Set rch_"+IntText(j+1)+" = [[
        Eval(\"Real observado_\"+label+\" = observado;\"),
        Eval(\"Real teorico_\"+label+\" = teorico;\")
      ]];")
    })
  };

  Serie uno = Not(residuos)+Not(Not(residuos));
  Serie not = Not(uno);

  Serie rch_pos = IfSer(GT(residuos,not),1,0);
  Serie rch_neg = IfSer(LT(residuos,not),1,0);

  Set sDtsRchPos = DatesRacha(rch_pos);
  Set sDtsRchNeg = DatesRacha(rch_neg);
  
  Set DistRacha(sDtsRchPos,residuos,"+")|DistRacha(sDtsRchNeg,residuos,"-")
};
//////////////////////////////////////////////////////////////////////////////
PutDescription("Devuelve el numero de rachas de residuos de igual signo 
consecutivos junto a la esperada teóricamente.",
Analisis.Rachas);
//////////////////////////////////////////////////////////////////////////////


//////////////////////////////////////////////////////////////////////////////
Set LookFor.VarianceChanges(Serie serie, Real long)
//////////////////////////////////////////////////////////////////////////////
{
//  Serie serie = var;
//  Real long   = 40;
  Date  iniDate = First(serie);
  Date  endDate = Last(serie);
  TimeSet tms   = Dating(serie);

  Real  card = Floor(
    SumS(Not(serie)+Not(Not(serie)))/long
  );

  Serie null = SubNullSer(tms,Succ(endDate,tms,-(long*card-1)),endDate);

  Set sPerform = For(1,card,Serie (Real j)
  {// Real j = 1;
    Date iniPerform = Succ(endDate,tms,-long*(j-1));
    Date endPerform = Succ(endDate,tms,-(long*j-1));
        
    //--> Las fechas están cambiadas bien.Es para empezar desde el final
    Serie sub = SubSer(serie,endPerform,iniPerform );
    Serie ci  = SubSer(CalInd(In(endPerform,iniPerform),tms),
      endPerform,
      iniPerform
    );

    //--> Saco la desviacion standard de los no nulos
    Real  sig = StDsS(AutoDating(sub));

    Serie end = sig*( null << ci >> null )
  });

  Serie sum = SetSumC(sPerform);
  Serie dif = (1-B):sum;

  //--> Elijo como candidatos a priori los mayores de dos sigmas
  Serie candidates = IfSer(LT(Serie Abs(dif),2*StDsS(AutoDating(dif))),0,dif);
  
  //--> La Kurtosis me dice si la distribucion es una normal, si es una 
  //    uniforme o en cambio es un pico 
  Real  rKurtosis = KurtosisS(AutoDating(dif));
  
  //--> La Kurtosis de una normal es 0, le doy 0.4 más para evitarme problemas
  //   de distribuciones muy parecidas a las normales
  Real rCriterium = 0 + 0.4;
  Set  sCandidates = If(Real Abs(rKurtosis)>rCriterium,
    DatesInd(candidates),
    Empty
  );

  //--> Me quedo con el último más long
  Date day = If( Card(sCandidates)>0,
    sCandidates[Card(sCandidates)],
    iniDate
  );

  Set  [[ sum,dif,candidates,day,rKurtosis ]]
};
//////////////////////////////////////////////////////////////////////////////
PutDescription("Busca cambios de varianza en una serie temporal. Recibe la 
serie que queramos estudiar y un numero real, que es un intervalo de dias. A
continuacion divide la serie en intervalos de longitud \"long\" y estudia los
cambios de varianza que se producen entre intervalos. Retorna un set con tres 
series. La primera muestra la varianza a intervalos (long), el segundo la 
serie diferenciada y la tercera la serie deiferenciada filtrada de candidatos 
pequeños las que informan sobre los posibles candidatos, junto con una fecha 
del dia más probable donde ocurra un cambio de varianza. Esto siempre y cuando
la distribucion de los candidatos sea más picuda que una normal, en tal caso se
considera que no hay cambio de varianza importante.",
LookFor.VarianceChanges);
//////////////////////////////////////////////////////////////////////////////


//////////////////////////////////////////////////////////////////////////////
Date VarianceChange(Serie serie, Real long)
//////////////////////////////////////////////////////////////////////////////
{
  Real iniTest = Floor(long/20);
  Real endTest = Floor(long/3);

  Set test = For(iniTest,endTest,Set (Real j)
  {
    Set research = LookFor.VarianceChanges(serie,10*j);
    Set proof    = If(Card(research)>0,research,Empty);

    Set PutName("Intervalo_"+IntText(j),proof)
  });

  Set sPulses = For(1,Card(test),Serie (Real k)
  {
    Date d = test[k][4];
    Serie SubSer(Pulse(d,Diario),First(serie),Last(serie))
  });

  Serie pulses = SetSum(sPulses);
  Serie select = IfSer(LT(Serie Abs(pulses),StDsS(AutoDating(pulses))),0,pulses);
  
  Set  sSelectcandidates = DatesInd(select);

  Date DayOfChange = If(Card(sSelectcandidates)>0,
    sSelectcandidates[Card(sSelectcandidates)],
    IniDate
  )
};
//////////////////////////////////////////////////////////////////////////////
PutDescription("Hace un estudio iterativo de la serie que le pasas, y devuelve
la fecha que estima que se ha producido un cambio de varianza. Si no detecta 
cambio alguno, devuelve la primera fecha de la serie. El campo \"long\" es el 
intervalo de dias en los que divide internamente la serie",
VarianceChange);
//////////////////////////////////////////////////////////////////////////////


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