root/library/bdm/estim/arx.h @ 649

/*!
  \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 "../math/functions.h"
#include "../stat/exp_family.h"
#include "../base/user_info.h"
//#include "../estim/kalman.h"
#include "arx_straux.h"

namespace bdm {

/*!
* \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]

See \ref tut_arx for mathematical treatment.

The easiest way how to use the class is:
\include arx_simple.cpp

        \todo sort out constant terms - bayes should accept vec without additional 1s
*/
class ARX: public BMEF {
protected:
        //!size of output variable (needed in regressors)
        int dimx;
        //!description of modelled data \f$ y_t \f$ in the likelihood function
        //! Do NOT access directly, only via \c get_yrv().
        RV _yrv;
        //! rv of regressor
        RV rgrrv;
        //! 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 ν
        //! switch if constant is modelled or not
        bool have_constant;
        //! cached value of data vector for have_constant =true
        vec _dt;
        //! Alternative estimate of parameters, used in stabilized forgetting, see [Kulhavy]
        egiw alter_est;
public:
        //! \name Constructors
        //!@{
        ARX ( const double frg0 = 1.0 ) : BMEF ( frg0 ), est (), V ( est._V() ), nu ( est._nu() ), have_constant(true), _dt() {};
        ARX ( const ARX &A0 ) : BMEF (A0.frg), est (A0.est), V ( est._V() ), nu ( est._nu() ), have_constant(A0.have_constant), _dt(A0._dt) {
                dimx = A0.dimx;
                _yrv = A0._yrv;
                rgrrv = A0.rgrrv;
                set_drv(A0._drv());
        };
        ARX* _copy_() const;
        void set_parameters ( double frg0 ) {
                frg = frg0;
        }
        void set_constant ( bool const0 ) {
                have_constant=const0;
        }
        void set_statistics ( int dimx0, const ldmat V0, double nu0 = -1.0 ) {
                est.set_parameters ( dimx0, V0, nu0 );
                last_lognc = est.lognc();
                dimx = dimx0;
        }
        //!@}

        //! Set sufficient statistics
        void set_statistics ( const BMEF* BM0 );

        //!\name Mathematical operations
        //!@{

        //! Weighted Bayes \f$ dt = [y_t psi_t] \f$.
        void bayes ( const vec &dt, const double w );
        void bayes ( const vec &dt ) {
                bayes ( dt, 1.0 );
        };
        double logpred ( const vec &dt ) const;
        void flatten ( const BMEF* B ) {
                const ARX* A = dynamic_cast<const ARX*> ( B );
                // nu should be equal to B.nu
                est.pow ( A->nu / nu );
                if ( evalll ) {
                        last_lognc = est.lognc();
                }
        }
        //! Conditioned version of the predictor
        enorm<ldmat>* epredictor ( const vec &rgr ) const;
        //! Predictor for empty regressor
        enorm<ldmat>* epredictor() const {
                bdm_assert_debug ( dimx == V.rows() - 1, "Regressor is not only 1" );
                return epredictor ( vec_1 ( 1.0 ) );
        }
        //! conditional version of the predictor
        mlnorm<ldmat>* predictor() const;
        mlstudent* predictor_student() const;
        //! Brute force structure estimation.\return indeces of accepted regressors.
        ivec structure_est ( egiw Eg0 );
        //! Smarter structure estimation by Ludvik Tesar.\return indeces of accepted regressors.
        ivec structure_est_LT ( egiw Eg0 );
        //!@}

        //!\name Access attributes
        //!@{
        const egiw& posterior() const {
                return est;
        }
        //!@}

        //!\name Connection
        //!@{
        void set_rv ( const RV &yrv0 , const RV &rgrrv0 ) {
                _yrv = yrv0;
                rgrrv=rgrrv0;
                set_drv(concat(yrv0, rgrrv));
        }

        RV& get_yrv() {
                //if yrv is not ready create it
                if ( _yrv._dsize() != dimx ) {
                        int i = 0;
                        while ( _yrv._dsize() < dimx ) {
                                _yrv.add ( drv ( vec_1 ( i ) ) );
                                i++;
                        }
                }
                //yrv should be ready by now
                bdm_assert_debug ( _yrv._dsize() == dimx, "incompatible drv" );
                return _yrv;
        }
        //!@}

        /*! UI for ARX estimator

        \code
        class = 'ARX';
        rv    = RV({names_of_dt} )                 // description of output variables
        rgr   = RV({names_of_regressors}, [-1,-2]} // description of regressor variables
        constant = 1;                              // 0/1 switch if the constant term is modelled or not

        --- optional ---
        V0  = [1 0;0 1];                           // Initial value of information matrix V
          --- OR ---
        dV0 = [1e-3, 1e-5, 1e-5, 1e-5];            // Initial value of diagonal of information matrix V
                                                                                           // default: 1e-3 for rv, 1e-5 for rgr
        nu0 = 6;                                                   // initial value of nu, default: rgrlen + 2
        frg = 1.0;                                 // forgetting, default frg=1.0

        rv_param   = RV({names_of_parameters}}     // description of parametetr names 
                                                                                           // default: ["theta_i" and "r_i"]
        \endcode
        */
        void from_setting ( const Setting &set );

        void validate() {
                bdm_assert(dimx == _yrv._dsize(), "RVs of parameters and regressor do not match");
                
        }
};

UIREGISTER ( ARX );
SHAREDPTR ( ARX );

/*! ARX model conditined by knowledge of the forgetting factor
\f[ f(\theta| d_1 \ldots d_t , \phi_t) \f]
*/
class ARXfrg : public ARX{
        public:
                ARXfrg():ARX(){};
                ARXfrg(const ARXfrg &A0):ARX(A0){};
                ARXfrg* _copy_() const {ARXfrg *A = new ARXfrg(*this); return A;}
        void condition(const vec &val){
                frg = val(0);
        }
};
UIREGISTER(ARXfrg);
};
#endif // AR_H