Changeset 369ee3d in tol


Ignore:
Timestamp:
Sep 26, 2012, 10:28:17 AM (14 years ago)
Author:
Víctor de Buen <vdebuen@…>
Branches:
1941-docker-image-to-test, devel, master
Children:
3cf31c1
Parents:
d404d85
Message:

Refs #1136
New methods to handle with conditioned least squares estimation and stationarity probability calculation

git-svn-id: ​https://www.tol-project.org/svn/tolp/trunk@5145 be4846ae-f051-483a-831f-b64e7d08527f

Location:
tol/bmath/bstat
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • tol/bmath/bstat/arma.cpp

    rd404d85 r369ee3d  
    3636#include <tol/tol_bratgra.h>
    3737#include <tol/tol_bar.h>
     38#include <tol/tol_bprdist.h>
    3839
    3940//--------------------------------------------------------------------
    … …  
    11541155//--------------------------------------------------------------------
    11551156{
     1157  BInt i,j;
     1158  if(&factor_ == &factor) { return; }
    11561159  if(!factor.Size())
    11571160  {
    … …  
    11591162    grFact_   = BArray<BARIMAFactor>(0);
    11601163    acumFact_ = BArray<BARIMAFactor>(0);
     1164    coef_.AllocBuffer(0);
    11611165    d_ = 0;
    11621166    p_ = 0;
    … …  
    11641168    return;
    11651169  }
    1166   BInt i,j;
    11671170  factor_.AllocBuffer(factor.Size());
    1168   for(i=j=0;i<factor.Size();i++)
     1171  for(i=j=n_=0;i<factor.Size();i++)
    11691172  {
    11701173    if(factor(i).HasValue())
    11711174    {
    11721175      factor_(j++) = factor(i);
     1176      n_+=(factor(i).ar_.Size()-1)+(factor(i).ma_.Size()-1);
    11731177    }
    11741178  }
    … …  
    11771181    factor_(j++) = factor(0);
    11781182  }
    1179 
    11801183  factor_.ReallocBuffer(j);
    11811184  factor_.Sort(ArrayARIMAFactorCmp);
     1185  coef_.ReallocBuffer(n_);
     1186  FactorUpdated();
     1187}
     1188
     1189//--------------------------------------------------------------------
     1190void BARIMA::FactorUpdated()
     1191//--------------------------------------------------------------------
     1192{
     1193  BInt i,j;
     1194  GetCoef();
    11821195  grFact_   = BArray<BARIMAFactor>(factor_.Size());
    11831196  acumFact_ = BArray<BARIMAFactor>(factor_.Size());
    … …  
    12121225  p_ = prod_.ar_ .Degree();
    12131226  q_ = prod_.ma_ .Degree();
    1214 //Std(BText("\nd=")+d_);
    1215 //Std(BText("\np=")+p_);
    1216 //Std(BText("\nq=")+q_);
    1217 
     1227}
     1228
     1229//--------------------------------------------------------------------
     1230const BArray<BDat>& BARIMA::GetCoef() const
     1231//--------------------------------------------------------------------
     1232{
     1233  int i,j,k,r;
     1234  for(r=k=0;k<factor_.Size();k++)
     1235  {
     1236    BPolyn<BDat>& ar = factor_(k).ar_;
     1237    BPolyn<BDat>& ma = factor_(k).ma_;
     1238    for(i=1;i<ar.Size();i++)
     1239    {
     1240      coef_[r++] = ar[i].Coef();
     1241    }
     1242    for(j=1;j<ma.Size();j++)
     1243    {
     1244      coef_[r++] = ma[j].Coef();
     1245    }
     1246  }
     1247  return(coef_);
     1248}
     1249
     1250//--------------------------------------------------------------------
     1251void BARIMA::PutCoef(const BArray<BDat>& coef)
     1252//--------------------------------------------------------------------
     1253{
     1254  int i,j,k,r;
     1255  if(&coef_ != &coef) { coef_ = coef; }
     1256  for(r=k=0;k<factor_.Size();k++)
     1257  {
     1258    BPolyn<BDat>& ar = factor_(k).ar_;
     1259    BPolyn<BDat>& ma = factor_(k).ma_;
     1260    for(i=1;i<ar.Size();i++)
     1261    {
     1262      ar[i].PutCoef(coef[r++]);
     1263    }
     1264    for(j=1;j<ma.Size();j++)
     1265    {
     1266      ma[j].PutCoef(coef[r++]);
     1267    }
     1268  }
     1269  FactorUpdated();
    12181270}
    12191271
    … …  
    12341286  BMat z1 = z_.Sub(d_,0,z_.Rows()-d_,z_.Columns());
    12351287  Backward(prod_.dif_,z0,z1,w_);
     1288  m_ = w_.Rows();
    12361289//TRZ(FrobeniusNormU(z0.Data())); TRZ(FrobeniusNormU(z1.Data()));TRZ(FrobeniusNormU(w_.Data()));
    12371290//TRZ(w_);
    … …  
    12441297}
    12451298
     1299//--------------------------------------------------------------------
     1300void BARIMA::CalcResiduals()
     1301//--------------------------------------------------------------------
     1302{
     1303  MatBackwardDifEq(prod_.ar_/prod_.ma_,w0_,w_,a0_,a_);
     1304}
     1305
     1306//--------------------------------------------------------------------
     1307void BARIMA::CalcResidualsJacobian(BMatrix<BDat>& J)
     1308//--------------------------------------------------------------------
     1309{
     1310  int i,r,k,j,d;
     1311  BMatrix<BDat> b;
     1312 
     1313  J.Alloc(m_,n_);
     1314  for(r=k=0;k<factor_.Size();k++)
     1315  {
     1316    BPolyn<BDat>& ar = factor_(k).ar_;
     1317    BPolyn<BDat>& ma = factor_(k).ma_;
     1318   
     1319    for(i=1;i<ar.Size();i++)
     1320    {
     1321      d = ar[i].Degree();
     1322      for(j=0; j<m_; j++) { J(j,k) = (j<d)?0:a_(j-d,0); }
     1323      k++;
     1324    }
     1325    MatBackwardDifEq(BPol::One()/ma,a_,b);
     1326    for(j=1;j<ma.Size();j++)
     1327    {
     1328      d = ma[i].Degree();
     1329      for(j=0; j<m_; j++) { J(j,k) = (j<d)?0:-b(j-d,0); }
     1330      k++;
     1331    }
     1332  }
     1333}
    12461334
    12471335//--------------------------------------------------------------------
    … …  
    12801368
    12811369//--------------------------------------------------------------------
    1282   bool BARIMA::CheckStationary()
     1370  bool BARIMA::CheckStationary(bool allowNonStat)
    12831371//--------------------------------------------------------------------
    12841372{
    … …  
    12941382  }
    12951383  bool ok = arOk && maOk;
    1296   if(!ok)
     1384  if(!allowNonStat && !ok)
    12971385  {
    12981386    Warning(I2("ARMA polynomials are not all stationary.",
    … …  
    13601448    if(calcInitValues && mpq)
    13611449    {
    1362       int j,d;
     1450      int j;
    13631451      BMatrix<BDat> covWW0(mpq,N);
    13641452      for(i=0; i<mpq; i++)
    … …  
    15861674}
    15871675
     1676
     1677
     1678//------------------------------------------------------------------
     1679BARIMACondLeastSqr::BARIMACondLeastSqr(BARIMA* model)
     1680: BRnRmFunction(), M_(model)
     1681//------------------------------------------------------------------
     1682{
     1683  SetDimensions(M_->n_, M_->m_);
     1684};
     1685
     1686//------------------------------------------------------------------
     1687void BARIMACondLeastSqr::Evaluate(      BArray<BDat>& a,
     1688              const BArray<BDat>& p)
     1689//------------------------------------------------------------------
     1690{
     1691  M_->PutCoef(p);
     1692  M_->CalcResiduals();
     1693  a = M_->a_.Data();
     1694};
     1695
     1696//------------------------------------------------------------------
     1697void BARIMACondLeastSqr::Jacobian(const BArray<BDat>&  x,
     1698              const BArray<BDat>&  y,
     1699                    BMatrix<BDat>& J)
     1700//------------------------------------------------------------------
     1701{
     1702  M_->PutCoef(x);
     1703  M_->CalcResidualsJacobian(J);
     1704};
     1705
     1706//------------------------------------------------------------------
     1707BBool BARIMACondLeastSqr::Marquardt()
     1708//------------------------------------------------------------------
     1709{
     1710  BDiagMatrix<BDat> Dp, D;
     1711  BMatrix<BDat> U,V,Vt;
     1712  x.Alloc(n_,1);
     1713  G.Alloc(n_,1);
     1714  H.Alloc(n_,n_);
     1715  U.Alloc(m_, n_);
     1716  V.Alloc(n_, n_);
     1717  D.Alloc(n_);
     1718  M_->stdErr_ = LeastSqrMarquardt(M_->coef_, M_->a_.GetData(), M_->Ja_);
     1719  stdErr = M_->stdErr_;
     1720  BDat var = stdErr^2;
     1721  x.GetData() = M_->coef_;
     1722  gsl_SingularValueDecomposition(M_->Ja_, U, D, V,
     1723  BText("Golub_Reinsch_Mod"));
     1724  Vt  = V.T();
     1725  Dp  = D.P(Sqrt(D.Rows()*BDat::Zero()));
     1726  H = (V*D*D*Vt)/var;
     1727  S = (V*Dp*Dp*Vt)*var;
     1728  SLt = V*Dp*stdErr;
     1729  return(stdErr.IsKnown());
     1730};
     1731
     1732//------------------------------------------------------------------
     1733BDat BARIMACondLeastSqr::StationarityProb(int sampleLength)
     1734//------------------------------------------------------------------
     1735{
     1736  int iter, t, nonStationaryCount=0;
     1737  BMatrix<BDat> e(n_,1), z(n_,1);
     1738  for(iter=1; iter<= sampleLength; iter++)
     1739  {
     1740    for(t=0; t<n_; t++)
     1741    {
     1742      e(t,0) = BNormalDist::Random01();
     1743    }
     1744    z = x + SLt * e;
     1745    M_->PutCoef(z.Data());
     1746    nonStationaryCount += !M_->CheckStationary(true); 
     1747  }
     1748  return(double(nonStationaryCount)/sampleLength);
     1749};
     1750
     1751
     1752/* */

  • tol/bmath/bstat/tol_barma.h

    rd404d85 r369ee3d  
    146146  BArray<BARIMAFactor>  grFact_;
    147147  BArray<BARIMAFactor>  acumFact_;
    148   BARIMAFactor          prod_;
     148  BARIMAFactor prod_;
     149  BArray<BDat> coef_;
    149150  BInt                  d_;
    150151  BInt                  p_;
    151152  BInt                  q_;
     153  BInt                  n_;
     154  BInt                  m_;
    152155  BMatrix<BDat>         z_;
    153156  BMatrix<BDat>         w_;
    … …  
    157160  BMatrix<BDat>         a0_;
    158161  BMatrix<BDat>         w0_;
     162  BMatrix<BDat>         Ja_;
    159163  BArray<BDat> aCov_;
    160164  BArray<BDat> aCor_;
    … …  
    167171  BDat likelihoodCoef_;
    168172  BDat waCoef_;
     173  BDat stdErr_;
    169174
    170175  BARIMA() {};
    171176 ~BARIMA() {};
    172177
     178  void PutOutputData(const BMatrix<BDat>& z);
    173179  void PutFactors   (const BArray<BARIMAFactor>& factor);
    174   void PutOutputData(const BMatrix<BDat>& z);
     180  void FactorUpdated();
     181  void PutCoef(const BArray<BDat>& coef);
     182  const BArray<BDat>& GetCoef() const;
    175183  void OutputDataUpdated();
    176184  void SetError();
    177   bool CheckStationary();
     185  void CalcResiduals();
     186  void CalcResidualsJacobian(BMatrix<BDat>& J);
     187  bool CheckStationary(bool allowNonStat=false);
    178188  bool CalcAutoCovarianze(BInt N);
    179189  bool CalcAutoCorrelation();
    … …  
    184194
    185195
     196//--------------------------------------------------------------------
     197class BARIMACondLeastSqr : public BRnRmFunction
     198
     199/*! Class to estimate an ARIMA model conditioned to initial values.
     200 */
     201//--------------------------------------------------------------------
     202{
     203public:
     204  BARIMA* M_;
     205  BMatrix<BDat> x;
     206  BMatrix<BDat> G;
     207  BMatrix<BDat> H;
     208  BMatrix<BDat> S;
     209  BMatrix<BDat> U;
     210  BDiagMatrix<BDat> D;
     211  BDiagMatrix<BDat> Dp;
     212  BMatrix<BDat> V;
     213  BMatrix<BDat> SLt;
     214  BMatrix<BDat> x_StdErr;
     215  BMatrix<BDat> x_T;
     216  BMatrix<BDat> x_RefuseProb;
     217  BDat stdErr;
     218  BDat S_logDet;
     219
     220  BARIMACondLeastSqr(BARIMA* model);
     221  void Evaluate(      BArray<BDat>& a,
     222                const BArray<BDat>& p);
     223  void Jacobian(const BArray<BDat>&  x,
     224                const BArray<BDat>&  y,
     225                      BMatrix<BDat>& J);
     226  BBool Marquardt();
     227  BDat StationarityProb(int sampleLength);
     228};
    186229#endif  //BARMA
Note: See TracChangeset for help on using the changeset viewer.