Index: bdm/estim/arx.h
===================================================================
--- bdm/estim/arx.h (revision 97)
+++ bdm/estim/arx.h (revision 97)
@@ -0,0 +1,61 @@
+/*!
+  \file
+  \brief Bayesian Filtering for generalized autoregressive (ARX) model
+  \author Vaclav Smidl.
+
+  -----------------------------------
+  BDM++ - C++ library for Bayesian Decision Making under Uncertainty
+
+  Using IT++ for numerical operations
+  -----------------------------------
+*/
+
+#ifndef AR_H
+#define AR_H
+
+#include <itpp/itbase.h>
+#include "../stat/libFN.h"
+#include "../stat/libEF.h"
+
+using namespace itpp;
+
+/*!
+* \brief Linear Autoregressive model with Gaussian noise
+
+Regression of the following kind:
+\f[
+y_t = \theta_1 \psi_1 + \theta_2 + \psi_2 +\ldots + \theta_n \psi_n + r e_t
+\f]
+where unknown parameters \c rv are \f$[\theta r]\f$, regression vector \f$\psi=\psi(y_{1:t},u_{1:t})\f$ is a known function of past outputs and exogeneous variables \f$u_t\f$. Distrubances \f$e_t\f$ are supposed to be normally distributed:
+\f[
+e_t \sim \mathcal{N}(0,1).
+\f]
+
+Extension for time-variant parameters \f$\theta_t,r_t\f$ may be achived using exponential forgetting (Kulhavy and Zarrop, 1993). In such a case, the forgetting factor \c frg \f$\in <0,1>\f$ should be given in the constructor. Time-invariant parameters are estimated for \c frg = 1.
+*/
+class ARX: public BM {
+protected:
+	//! Posterior estimate of \f$\theta,r\f$ in the form of Normal-inverse Wishart density
+	egiw est;
+	//! cached value of est.V
+	ldmat &V;
+	//! cached value of est.nu
+	double &nu;
+	//! forgetting factor 
+	double frg;
+	//! cached value of lognc() in the previous step (used in evaluation of \c ll )
+	double last_lognc;
+public:
+	//! Full constructor
+	ARX (RV &rv, mat &V0, double &nu0, double frg0=1.0) : BM(rv),est(rv,V0,nu0), V(est._V()), nu(est._nu()), frg(frg0){last_lognc=est.lognc();};
+	//! Here \f$dt = [y_t psi_t] \f$.
+	void bayes ( const vec &dt );
+	epdf& _epdf() {return est;}
+	//! Brute force structure estimation.\return indeces of accepted regressors.
+	ivec structure_est(egiw Eg0);
+};
+
+
+#endif // AR_H
+
+
Index: bdm/estim/arx.cpp
===================================================================
--- bdm/estim/arx.cpp (revision 14)
+++ bdm/estim/arx.cpp (revision 97)
@@ -1,6 +1,72 @@
-#include <itpp/itbase.h>
-#include "libBM.h"
+#include "arx.h"
 
 using namespace itpp;
 
+void ARX::bayes ( const vec &dt ) {
+	double lnc;
 
+	V.opupdt ( dt,frg );
+	nu+=frg;
+
+	if ( evalll ) {
+		lnc = est.lognc();
+		ll = lnc - last_lognc;
+		last_lognc = lnc;
+	}
+}
+
+/*! \brief Return the best structure
+@param Eg a copy of GiW density that is being examined
+@param Eg0 a copy of prior GiW density before estimation
+@param Egll likelihood of the current Eg
+@param indeces current indeces
+\return best likelihood in the structure below the given one
+*/
+double egiw_bestbelow ( egiw Eg, egiw Eg0, double Egll, ivec &indeces ) { //parameter Eg is a copy!
+	ldmat Vo = Eg._V(); //copy
+	ldmat Vo0 = Eg._V(); //copy
+	ldmat& Vp = Eg._V(); // pointer into Eg
+	ldmat& Vp0 = Eg._V(); // pointer into Eg
+	int end = Vp.rows()-1;
+	int i;
+	mat Li;
+	mat Li0;
+	double maxll=Egll;
+	double tmpll=Egll;
+	double belll=Egll;
+
+	ivec tmpindeces;
+	ivec maxindeces=indeces;
+
+	//try to remove only one rv
+	for ( i=0;i<end;i++ ) {
+		//copy original
+		Li = Vo._L();
+		Li0 = Vo0._L();
+		//remove stuff
+		Li.del_col ( i+1 );
+		Li0.del_col ( i+1 );
+		Vp.ldform ( Li,Vo._D() );
+		Vp0.ldform ( Li0,Vo0._D() );
+		tmpll = Eg.lognc()-Eg0.lognc(); // likelihood is difference of norm. coefs.
+		//
+		if ( tmpll > Egll ) { //increase of the likelihood
+			tmpindeces = indeces;
+			tmpindeces.del ( i );
+			//search for a better match in this substructure
+			belll=egiw_bestbelow ( Eg, Eg0, tmpll, tmpindeces );
+			if ( belll>maxll ) { //better match found
+				maxll = belll;
+				maxindeces = tmpindeces;
+			}
+		}
+	}
+	indeces = maxindeces;
+	return maxll;
+}
+
+ivec ARX::structure_est ( egiw est0 ) {
+	ivec ind=linspace ( 1,rv.count()-1 );
+	double tmp = egiw_bestbelow ( est, est0, est.lognc()- est0.lognc(), ind );
+	return ind;
+}
