Index: tol/bmath/bstat/arma.cpp
===================================================================
--- tol/bmath/bstat/arma.cpp	(revision d404d8572fbdee492d724b1d45822205c08e2998)
+++ tol/bmath/bstat/arma.cpp	(revision 369ee3d25c07f1908533cacf8a2a8830cd37cc50)
@@ -36,4 +36,5 @@
 #include <tol/tol_bratgra.h>
 #include <tol/tol_bar.h>
+#include <tol/tol_bprdist.h>
 
 //--------------------------------------------------------------------
@@ -1154,4 +1155,6 @@
 //--------------------------------------------------------------------
 {
+  BInt i,j;
+  if(&factor_ == &factor) { return; }
   if(!factor.Size())
   {
@@ -1159,4 +1162,5 @@
     grFact_   = BArray<BARIMAFactor>(0);
     acumFact_ = BArray<BARIMAFactor>(0);
+    coef_.AllocBuffer(0);
     d_ = 0;
     p_ = 0;
@@ -1164,11 +1168,11 @@
     return;
   }
-  BInt i,j;
   factor_.AllocBuffer(factor.Size());
-  for(i=j=0;i<factor.Size();i++)
+  for(i=j=n_=0;i<factor.Size();i++)
   {
     if(factor(i).HasValue())
     {
       factor_(j++) = factor(i);
+      n_+=(factor(i).ar_.Size()-1)+(factor(i).ma_.Size()-1);
     }
   }
@@ -1177,7 +1181,16 @@
     factor_(j++) = factor(0);
   }
-
   factor_.ReallocBuffer(j);
   factor_.Sort(ArrayARIMAFactorCmp);
+  coef_.ReallocBuffer(n_);
+  FactorUpdated();
+}
+
+//--------------------------------------------------------------------
+void BARIMA::FactorUpdated()
+//--------------------------------------------------------------------
+{
+  BInt i,j;
+  GetCoef();
   grFact_   = BArray<BARIMAFactor>(factor_.Size());
   acumFact_ = BArray<BARIMAFactor>(factor_.Size());
@@ -1212,8 +1225,47 @@
   p_ = prod_.ar_ .Degree();
   q_ = prod_.ma_ .Degree();
-//Std(BText("\nd=")+d_);
-//Std(BText("\np=")+p_);
-//Std(BText("\nq=")+q_);
-
+}
+
+//--------------------------------------------------------------------
+const BArray<BDat>& BARIMA::GetCoef() const
+//--------------------------------------------------------------------
+{
+  int i,j,k,r;
+  for(r=k=0;k<factor_.Size();k++)
+  {
+    BPolyn<BDat>& ar = factor_(k).ar_;
+    BPolyn<BDat>& ma = factor_(k).ma_;
+    for(i=1;i<ar.Size();i++)
+    {
+      coef_[r++] = ar[i].Coef();
+    }
+    for(j=1;j<ma.Size();j++)
+    {
+      coef_[r++] = ma[j].Coef();
+    }
+  }
+  return(coef_);
+}
+
+//--------------------------------------------------------------------
+void BARIMA::PutCoef(const BArray<BDat>& coef)
+//--------------------------------------------------------------------
+{
+  int i,j,k,r;
+  if(&coef_ != &coef) { coef_ = coef; }
+  for(r=k=0;k<factor_.Size();k++)
+  {
+    BPolyn<BDat>& ar = factor_(k).ar_;
+    BPolyn<BDat>& ma = factor_(k).ma_;
+    for(i=1;i<ar.Size();i++)
+    {
+      ar[i].PutCoef(coef[r++]);
+    }
+    for(j=1;j<ma.Size();j++)
+    {
+      ma[j].PutCoef(coef[r++]);
+    }
+  }
+  FactorUpdated();
 }
 
@@ -1234,4 +1286,5 @@
   BMat z1 = z_.Sub(d_,0,z_.Rows()-d_,z_.Columns());
   Backward(prod_.dif_,z0,z1,w_);
+  m_ = w_.Rows();
 //TRZ(FrobeniusNormU(z0.Data())); TRZ(FrobeniusNormU(z1.Data()));TRZ(FrobeniusNormU(w_.Data()));
 //TRZ(w_);
@@ -1244,4 +1297,39 @@
 }
 
+//--------------------------------------------------------------------
+void BARIMA::CalcResiduals()
+//--------------------------------------------------------------------
+{
+  MatBackwardDifEq(prod_.ar_/prod_.ma_,w0_,w_,a0_,a_);
+}
+
+//--------------------------------------------------------------------
+void BARIMA::CalcResidualsJacobian(BMatrix<BDat>& J)
+//--------------------------------------------------------------------
+{
+  int i,r,k,j,d;
+  BMatrix<BDat> b;
+  
+  J.Alloc(m_,n_);
+  for(r=k=0;k<factor_.Size();k++)
+  {
+    BPolyn<BDat>& ar = factor_(k).ar_;
+    BPolyn<BDat>& ma = factor_(k).ma_;
+    
+    for(i=1;i<ar.Size();i++)
+    {
+      d = ar[i].Degree();
+      for(j=0; j<m_; j++) { J(j,k) = (j<d)?0:a_(j-d,0); } 
+      k++;
+    } 
+    MatBackwardDifEq(BPol::One()/ma,a_,b);
+    for(j=1;j<ma.Size();j++)
+    {
+      d = ma[i].Degree();
+      for(j=0; j<m_; j++) { J(j,k) = (j<d)?0:-b(j-d,0); } 
+      k++;
+    }
+  }
+}
 
 //--------------------------------------------------------------------
@@ -1280,5 +1368,5 @@
 
 //--------------------------------------------------------------------
-  bool BARIMA::CheckStationary()
+  bool BARIMA::CheckStationary(bool allowNonStat)
 //--------------------------------------------------------------------
 {
@@ -1294,5 +1382,5 @@
   }
   bool ok = arOk && maOk;
-  if(!ok)
+  if(!allowNonStat && !ok)
   {
     Warning(I2("ARMA polynomials are not all stationary.",
@@ -1360,5 +1448,5 @@
     if(calcInitValues && mpq)
     { 
-      int j,d;
+      int j;
       BMatrix<BDat> covWW0(mpq,N);
       for(i=0; i<mpq; i++) 
@@ -1586,2 +1674,79 @@
 }
 
+
+
+//------------------------------------------------------------------
+BARIMACondLeastSqr::BARIMACondLeastSqr(BARIMA* model)
+: BRnRmFunction(), M_(model)
+//------------------------------------------------------------------
+{
+  SetDimensions(M_->n_, M_->m_);
+};
+
+//------------------------------------------------------------------
+void BARIMACondLeastSqr::Evaluate(      BArray<BDat>& a,
+              const BArray<BDat>& p)
+//------------------------------------------------------------------
+{
+  M_->PutCoef(p);
+  M_->CalcResiduals();
+  a = M_->a_.Data();
+};
+
+//------------------------------------------------------------------
+void BARIMACondLeastSqr::Jacobian(const BArray<BDat>&  x,
+              const BArray<BDat>&  y,
+                    BMatrix<BDat>& J)
+//------------------------------------------------------------------
+{
+  M_->PutCoef(x);
+  M_->CalcResidualsJacobian(J);
+};
+
+//------------------------------------------------------------------
+BBool BARIMACondLeastSqr::Marquardt()
+//------------------------------------------------------------------
+{
+  BDiagMatrix<BDat> Dp, D;
+  BMatrix<BDat> U,V,Vt;
+  x.Alloc(n_,1);
+  G.Alloc(n_,1);
+  H.Alloc(n_,n_);
+  U.Alloc(m_, n_);
+  V.Alloc(n_, n_);
+  D.Alloc(n_);
+  M_->stdErr_ = LeastSqrMarquardt(M_->coef_, M_->a_.GetData(), M_->Ja_);
+  stdErr = M_->stdErr_;
+  BDat var = stdErr^2;
+  x.GetData() = M_->coef_;
+  gsl_SingularValueDecomposition(M_->Ja_, U, D, V,
+  BText("Golub_Reinsch_Mod"));
+  Vt  = V.T();
+  Dp  = D.P(Sqrt(D.Rows()*BDat::Zero()));
+  H = (V*D*D*Vt)/var;
+  S = (V*Dp*Dp*Vt)*var;
+  SLt = V*Dp*stdErr;
+  return(stdErr.IsKnown());
+};
+
+//------------------------------------------------------------------
+BDat BARIMACondLeastSqr::StationarityProb(int sampleLength)
+//------------------------------------------------------------------
+{
+  int iter, t, nonStationaryCount=0;
+  BMatrix<BDat> e(n_,1), z(n_,1);
+  for(iter=1; iter<= sampleLength; iter++)
+  {
+    for(t=0; t<n_; t++)
+    {
+      e(t,0) = BNormalDist::Random01();
+    }
+    z = x + SLt * e;
+    M_->PutCoef(z.Data());
+    nonStationaryCount += !M_->CheckStationary(true);  
+  }
+  return(double(nonStationaryCount)/sampleLength);
+};
+
+
+/* */
Index: tol/bmath/bstat/tol_barma.h
===================================================================
--- tol/bmath/bstat/tol_barma.h	(revision d404d8572fbdee492d724b1d45822205c08e2998)
+++ tol/bmath/bstat/tol_barma.h	(revision 369ee3d25c07f1908533cacf8a2a8830cd37cc50)
@@ -146,8 +146,11 @@
   BArray<BARIMAFactor>	grFact_;
   BArray<BARIMAFactor>	acumFact_;
-  BARIMAFactor		prod_;
+  BARIMAFactor prod_;
+  BArray<BDat> coef_;
   BInt			d_;
   BInt			p_;
   BInt			q_;
+  BInt			n_;
+  BInt			m_;
   BMatrix<BDat>		z_;
   BMatrix<BDat>		w_;
@@ -157,4 +160,5 @@
   BMatrix<BDat>		a0_;
   BMatrix<BDat>		w0_;
+  BMatrix<BDat>		Ja_;
   BArray<BDat> aCov_;
   BArray<BDat> aCor_;
@@ -167,13 +171,19 @@
   BDat likelihoodCoef_;
   BDat waCoef_;
+  BDat stdErr_;
 
   BARIMA() {};
  ~BARIMA() {};
 
+  void PutOutputData(const BMatrix<BDat>& z);
   void PutFactors   (const BArray<BARIMAFactor>& factor);
-  void PutOutputData(const BMatrix<BDat>& z);
+  void FactorUpdated();
+  void PutCoef(const BArray<BDat>& coef);
+  const BArray<BDat>& GetCoef() const;
   void OutputDataUpdated();
   void SetError();
-  bool CheckStationary();
+  void CalcResiduals();
+  void CalcResidualsJacobian(BMatrix<BDat>& J);
+  bool CheckStationary(bool allowNonStat=false);
   bool CalcAutoCovarianze(BInt N);
   bool CalcAutoCorrelation();
@@ -184,3 +194,36 @@
 
 
+//--------------------------------------------------------------------
+class BARIMACondLeastSqr : public BRnRmFunction
+
+/*! Class to estimate an ARIMA model conditioned to initial values.
+ */
+//--------------------------------------------------------------------
+{
+public:
+  BARIMA* M_;
+  BMatrix<BDat> x;
+  BMatrix<BDat> G;
+  BMatrix<BDat> H;
+  BMatrix<BDat> S;
+  BMatrix<BDat> U;
+  BDiagMatrix<BDat> D;
+  BDiagMatrix<BDat> Dp;
+  BMatrix<BDat> V;
+  BMatrix<BDat> SLt;
+  BMatrix<BDat> x_StdErr;
+  BMatrix<BDat> x_T;
+  BMatrix<BDat> x_RefuseProb;
+  BDat stdErr;
+  BDat S_logDet;
+
+  BARIMACondLeastSqr(BARIMA* model);
+  void Evaluate(      BArray<BDat>& a,
+                const BArray<BDat>& p);
+  void Jacobian(const BArray<BDat>&  x,
+                const BArray<BDat>&  y,
+                      BMatrix<BDat>& J);
+  BBool Marquardt();
+  BDat StationarityProb(int sampleLength);
+};
 #endif	//BARMA