Index: tol/bmath/bstat/arma.cpp
===================================================================
--- tol/bmath/bstat/arma.cpp	(revision 516e91c7aaba3a972d8d24b66ef89fa3013d6021)
+++ tol/bmath/bstat/arma.cpp	(revision 7168a546902fe27f45303698b15b71704584b420)
@@ -25,4 +25,5 @@
 
 #include <tol/tol_matrix.h>
+#include <tol/tol_bvmat.h>
 #include <tol/tol_barma.h>
 #include <tol/tol_bstat.h>
@@ -579,8 +580,17 @@
     BMat gp, test;
     BDat err;
+    
+    BVMat ms, As, Ls, gps;
+    As.DMat2sparse((const BMatrix<double>&)A);
+    ms.DMat2sparse((const BMatrix<double>&)m);
+    BVMat::CholeskiFactor(As,Ls,"XtX",true,false);
+    BVMat::CholeskiSolve(Ls,(ms.T()*As).T(),gps,"PtLLtP");
+    gps.GetDMat((BMatrix<double>&)gp);
+    test = A*gp-m;
+    err  = test.FrobeniusNorm();
 
 //*
   //TRZ(psi); TRZ(m); TRZ(A);
-
+/*
     if(p>100)
     {
@@ -602,6 +612,9 @@
       }
     }
+*/
+    bool closeToNonStationary = false;
     if(err.IsUnknown()  || (err>Sqrt((1+p)*DEpsilon()))) 
     {
+      closeToNonStationary = true;
       BText method = BText("Golub_Reinsch_Mod");
       BMat U, V;
@@ -622,6 +635,6 @@
       return(false);
     }
-    else if(err>Sqrt((1+p)*DEpsilon()))
-    {  
+    else if(closeToNonStationary)
+    { 
       Warning(I2("ARMA structure is not numerically stationary. "
                  "Lost precision for autocovarianze matrix: ",