//+------------------------------------------------------------------+
//|                                                RLSRegression.mqh |
//+------------------------------------------------------------------+
#ifndef RLSREGRESSION_MQH
#define RLSREGRESSION_MQH

//+------------------------------------------------------------------+
//| Class CRLSRegression                                             |
//| Purpose: maintains a recursively updated linear regression of    |
//| y on a time index, using the Sherman-Morrison rank-1 update      |
//| with a forgetting factor for non-stationary tracking.            |
//+------------------------------------------------------------------+
class CRLSRegression
  {
private:
   double            m_theta0;
   double            m_theta1;
   double            m_P00;
   double            m_P01;
   double            m_P10;
   double            m_P11;
   double            m_lambda;
   double            m_delta;
   int               m_obs_count;
   int               m_min_obs;

public:
                     CRLSRegression(void);
                    ~CRLSRegression(void);

   void              Init(double lambda,double delta,int min_obs);
   void              Reset(void);
   void              Update(double y,double t_index);
   double            GetSlope(void) const;
   double            GetIntercept(void) const;
   double            Forecast(double t_index) const;
   bool              IsValid(void) const;
   int               GetObsCount(void) const;
  };

//+------------------------------------------------------------------+
//| Constructor                                                      |
//+------------------------------------------------------------------+
CRLSRegression::CRLSRegression(void)
  {
   m_lambda=0.98;
   m_delta=1000.0;
   m_min_obs=20;
   m_obs_count=0;
   m_theta0=0.0;
   m_theta1=0.0;
   m_P00=m_delta;
   m_P01=0.0;
   m_P10=0.0;
   m_P11=m_delta;
  }

//+------------------------------------------------------------------+
//| Destructor                                                       |
//+------------------------------------------------------------------+
CRLSRegression::~CRLSRegression(void)
  {
  }

//+------------------------------------------------------------------+
//| Configure the filter's parameters and reset its state            |
//+------------------------------------------------------------------+
void CRLSRegression::Init(double lambda,double delta,int min_obs)
  {
   m_lambda=lambda;
   m_delta=delta;
   m_min_obs=min_obs;
   Reset();
  }

//+------------------------------------------------------------------+
//| Reset theta to zero and P to delta times the identity matrix     |
//+------------------------------------------------------------------+
void CRLSRegression::Reset(void)
  {
   m_theta0=0.0;
   m_theta1=0.0;
   m_P00=m_delta;
   m_P01=0.0;
   m_P10=0.0;
   m_P11=m_delta;
   m_obs_count=0;
  }

//+------------------------------------------------------------------+
//| Apply the Sherman-Morrison rank-1 RLS update for one new         |
//| observation y at time index t_index                              |
//+------------------------------------------------------------------+
void CRLSRegression::Update(double y,double t_index)
  {
//--- feature vector is x = [1, t_index]
   double x0=1.0;
   double x1=t_index;

//--- P_{t-1} * x_t  (a 2x1 vector, computed from the 2x2 P and x)
   double Px0=m_P00*x0+m_P01*x1;
   double Px1=m_P10*x0+m_P11*x1;

//--- denominator: lambda + x_t' * P_{t-1} * x_t
   double denom=m_lambda+(x0*Px0+x1*Px1);

//--- gain vector K_t = (P_{t-1} * x_t) / denom
   double K0=Px0/denom;
   double K1=Px1/denom;

//--- innovation: prediction error using the previous theta
   double y_hat=m_theta0*x0+m_theta1*x1;
   double innovation=y-y_hat;

//--- coefficient update: theta_t = theta_{t-1} + K_t * innovation
   m_theta0=m_theta0+K0*innovation;
   m_theta1=m_theta1+K1*innovation;

//--- x_t' * P_{t-1}  (a 1x2 row vector)
   double xP0=x0*m_P00+x1*m_P10;
   double xP1=x0*m_P01+x1*m_P11;

//--- P update: P_t = (P_{t-1} - K_t * (x_t' * P_{t-1})) / lambda
   double newP00=(m_P00-K0*xP0)/m_lambda;
   double newP01=(m_P01-K0*xP1)/m_lambda;
   double newP10=(m_P10-K1*xP0)/m_lambda;
   double newP11=(m_P11-K1*xP1)/m_lambda;

   m_P00=newP00;
   m_P01=newP01;
   m_P10=newP10;
   m_P11=newP11;

   m_obs_count++;
  }

//+------------------------------------------------------------------+
//| Accessor: current slope estimate (theta1)                        |
//+------------------------------------------------------------------+
double CRLSRegression::GetSlope(void) const
  {
   return(m_theta1);
  }

//+------------------------------------------------------------------+
//| Accessor: current intercept estimate (theta0)                    |
//+------------------------------------------------------------------+
double CRLSRegression::GetIntercept(void) const
  {
   return(m_theta0);
  }

//+------------------------------------------------------------------+
//| Compute the model's forecast at a given time index               |
//+------------------------------------------------------------------+
double CRLSRegression::Forecast(double t_index) const
  {
   return(m_theta0+m_theta1*t_index);
  }

//+------------------------------------------------------------------+
//| True once at least m_min_obs observations have been processed    |
//+------------------------------------------------------------------+
bool CRLSRegression::IsValid(void) const
  {
   return(m_obs_count>=m_min_obs);
  }

//+------------------------------------------------------------------+
//| Accessor: number of observations processed so far                |
//+------------------------------------------------------------------+
int CRLSRegression::GetObsCount(void) const
  {
   return(m_obs_count);
  }

#endif // RLSREGRESSION_MQH
//+------------------------------------------------------------------+