GRFKalmanAG

#property copyright "Copyleft © 2007, GammaRat Forex"
#property link      "http://www.gammarat.com/Forex/"

   

#include <GRFMatrixMath.mqh>
#property indicator_chart_window
#property indicator_buffers 7
#property indicator_color1 Yellow
#property indicator_color2 Yellow
#property indicator_color3 Yellow
#property indicator_color4 Red
#property indicator_color5 Red
#property indicator_color6 Red
#property indicator_color7 Red
#property indicator_style1 0
#property indicator_style2 2
#property indicator_style3 2
#property indicator_style4 2
#property indicator_style5 2
#property indicator_style6 0
#property indicator_style7 0


//---- input parameters
extern double Samples=60;
extern double DevLevel1 = 1;
extern double DevLevel2 = 2;
extern double Suppression_dB = 0;
extern double PredictBars1 = 4;
extern double PredictBars2 = 8;
extern bool PrintCurrentState = true;
extern bool AutoGainOn = true;


//---- buffers
double KalmanBuffer[];
double KalmanBufferPlus1[];
double KalmanBufferNeg1[];
double KalmanBufferPlus2[];
double KalmanBufferNeg2[];
double KalmanBufferPlus3[];
double KalmanBufferNeg3[];
double xkk[2][1],xkk1[2][1],xk1k1[2][1],yk[1][1],zk[1][1];
double Pkk[2][2],Pkk1[2][2],Pk1k1[2][2],Pkk_inv[2][2];
double Qkk[2][2],Hk[1][2],Hk_t[2][1],Sk[1][1],Sk_inv[1][1],Rk[1][1],Kk[2][1];
double Fk[2][2],Fk_t[2][2],GGT[2][2];
double Eye[2][2];
double temp22[2][2],temp21[2][1],temp12[1][2],temp11[1][1];
double xp[1][1],xpt[1][1],xpt_t[1][1];
double Phi1[1][2],Phi2[1][2];
double xdel[1][1],xdel_t[1][1];
double LookAhead=0;
static bool initiated=false;
static double c0_above=0, c0_below = 0;
static double c1_above=0, c1_below = 0;
static int last_time=-1;
static int tick_count=0;
//+------------------------------------------------------------------+
//| Custom indicator initialization function                         |
//+------------------------------------------------------------------+
int init()
  {
//---- indicators
   if(LookAhead <0)LookAhead=0;
   SetIndexStyle(0,DRAW_LINE);
   SetIndexBuffer(0,KalmanBuffer);
   SetIndexShift(0,LookAhead);
   SetIndexDrawBegin(0,LookAhead);
   SetIndexLabel(0,"Kalman Trend");
   if(MathAbs(DevLevel1) > 0) {
      SetIndexStyle(1,DRAW_LINE);
      SetIndexBuffer(1,KalmanBufferPlus1);
      SetIndexShift(1,0);
      SetIndexDrawBegin(1,1);
      SetIndexLabel(1,"Kalman +" + DoubleToStr(DevLevel1,1)  + " STD");
      SetIndexStyle(2,DRAW_LINE);
      SetIndexBuffer(2,KalmanBufferNeg1);
      SetIndexShift(2,0);
      SetIndexDrawBegin(2,1);
      SetIndexLabel(2,"Kalman -" + DoubleToStr(DevLevel1,1) + " STD");
    }
   if(MathAbs(DevLevel2) > 0){
      SetIndexStyle(3,DRAW_LINE);
      SetIndexBuffer(3,KalmanBufferPlus2);
      SetIndexShift(3,0);
      SetIndexDrawBegin(3,1);
      SetIndexLabel(3,"Kalman +" + DoubleToStr(DevLevel2,1) + " STD");
      SetIndexStyle(4,DRAW_LINE);
      SetIndexBuffer(4,KalmanBufferNeg2);
      SetIndexShift(4,0);
      SetIndexDrawBegin(4,1);
      SetIndexLabel(4,"Kalman -" + DoubleToStr(DevLevel2,1) + " STD");
      
      SetIndexStyle(5,DRAW_LINE);
      SetIndexBuffer(5,KalmanBufferPlus3);
      SetIndexShift(5,0);
      SetIndexDrawBegin(5,1);
      SetIndexLabel(5,"Kalman High +" + DoubleToStr(DevLevel2,1) + " STD");
      SetIndexStyle(6,DRAW_LINE);
      SetIndexBuffer(6,KalmanBufferNeg3);
      SetIndexShift(6,0);
      SetIndexDrawBegin(6,1);
      SetIndexLabel(6,"Kalman Low -" + DoubleToStr(DevLevel2,1) + " STD");
   }
//----
   return(0);
  }
//| Point and figure                                |
//+------------------------------------------------------------------+
int start()
{
   compute();
   return(0);
}
int compute(){
   int    i,j,counted_bars=IndicatorCounted();
   
   double g,g1;
   double c0;
 
   static double acc_sigma2;
   static double z_sigma2;
   static double normalization_time;
   static double working_samples;
   double z,diff; 
   double ts[3][1];
   static double dev_level1,dev_level2; 
   double temp11a[1][1],temp21a[2][1];
   double temp22a[2][2]; 
   static double ZGain = 0;
   static double AGain=1;
   static double tgain=1;
   static double time_last;
   static double zk_saved[2][1];
   int trial;
   
   //----
   //give it some time to warm up
   if(Bars<Samples*2) 
      return(0);
   //rewrite the last bar
   //if(counted_bars > Bars-2)return(0);
   //Print(Time[0]);
   if(Time[0] <= last_time +1) return(0);
   tick_count--;
   if(tick_count>0)return(0);
   tick_count=2;
   last_time=Time[0];
   
   //if(counted_bars>0)counted_bars--;
   
   //if(counted_bars >2) counted_bars -= 2;
   
   if( !initiated){
       
      Phi1[0][0] = 1;
      Phi1[0][1] = PredictBars1;
      Phi2[0][0] = 1;
      Phi2[0][1] = PredictBars2;
            
      dev_level1 = DevLevel1*MathSqrt(2);
      dev_level2 = DevLevel2*MathSqrt(2);
      ArrayCopy(KalmanBuffer,High);
      xkk1[0][0] = get_avg(Bars-1);
      xkk1[1][0] = 0;
      //MatrixPrint(xkk);
      MatrixEye(Pkk1);
      Fk[0][0] = 1;
      Fk[0][1]=1;
      Fk[1][0]=0;
      Fk[1][1] =1;
      MatrixTranspose(Fk,Fk_t);
      Qkk[0][0] = 1;
      Qkk[0][1] = 0;
      Qkk[1][0]=0;
      Qkk[1][1] = 1;
      GGT[0][0]=.25;
      GGT[1][0]=.5;
      GGT[0][1]=.5;
      GGT[1][1]=1;
      Hk[0][0] = 1;
      Hk[0][1] = 0;
      MatrixTranspose(Hk,Hk_t);
      
      Rk[0][0] = .1;
      MatrixEye(Eye);
      ZGain = 0;
      z_sigma2 = MathPow(10,-ZGain);
      acc_sigma2 = Point*MathPow(10,-Suppression_dB/20.);//*MathSqrt(0.0001/Point);
      MatrixScalarMul(acc_sigma2,Pkk1);
      working_samples = Samples;
      acc_sigma2 *= acc_sigma2;
       Rk[0][0] = z_sigma2;
      initiated = true;
      MatrixScalarMul(acc_sigma2,Pkk);
      //Print(1./MathSqrt(acc_sigma2));
      //Print(" ");
      AGain=1;
      zk_saved[0][0] = get_avg(Bars-1);
      zk_saved[1][0] = 0;
      //double A1[20][2],A1t[2][20],A11[2][2],A11_inv[2][2],X[2],B[20];
      //for(i=0;i<20;i++){
         
   }
 
   for(i=Bars-counted_bars; i> 0;i--){ //ticks are handled elsewhere
     
     //Start the update procedure
      MatrixZero(yk);
      zk[0][0] = get_avg(i);
      zk[1][0]= zk[0][0]-zk_saved[0][0];
      ArrayCopy(zk_saved,zk);
      diff = (zk[0][0]*Point-KalmanBuffer[i]);
      if(diff>=0) {
         diff = diff*diff;
         c0_above = (c0_above*(working_samples-1)+MathPow(get_avg(i)*Point-KalmanBuffer[i],2))/working_samples;
         c0_below = c0_below*(working_samples-1)/working_samples;
         c1_above = (c1_above*(working_samples-1)+MathPow(High[i]-KalmanBuffer[i],2))/working_samples;
         c1_below = c1_below*(working_samples-1)/working_samples;
      }else{
         diff = diff*diff;
         c0_below = (c0_below*(working_samples-1)+MathPow(get_avg(i)*Point-KalmanBuffer[i],2))/working_samples;
         c0_above = c0_above*(working_samples-1)/working_samples;
         c1_below = (c1_below*(working_samples-1)+MathPow(Low[i]-KalmanBuffer[i],2))/working_samples;
         c1_above = c1_above*(working_samples-1)/working_samples;
      }  
      if(MathAbs(DevLevel1)>0){
         KalmanBufferPlus1[i-1] = KalmanBuffer[i]+dev_level1*MathSqrt(c0_above);
         KalmanBufferNeg1[i-1] = KalmanBuffer[i]-dev_level1*MathSqrt(c0_below);
      }
      if(MathAbs(DevLevel2)>0){   
         KalmanBufferPlus2[i-1] = KalmanBuffer[i]+dev_level2*MathSqrt(c0_above);
         KalmanBufferNeg2[i-1] = KalmanBuffer[i]-dev_level2*MathSqrt(c0_below);
         KalmanBufferPlus3[i-1] = KalmanBuffer[i]+dev_level2*MathSqrt(c1_above);
         KalmanBufferNeg3[i-1] = KalmanBuffer[i]-dev_level2*MathSqrt(c1_below);
      } 
      //this is always left alone, for now
           
     
      ArrayCopy(Qkk,GGT);
      
      // set the auto gain
      if(AutoGainOn && MathSqrt((c0_above+c0_below)/2)> Point){   
         if(c0_below >0 || c0_above>0)
            tgain = MathAbs((c0_below)/(c0_below+c0_above));
         else tgain = 0.5;
         if(tgain< 0.2 || tgain > 0.8)
            AGain = AGain*MathPow(10,0.005);
         else if(tgain> 0.4 && tgain < 0.6)
            AGain=AGain*MathPow(10,-0.005);
      }else{
         AGain = 1;
      }   
      tgain = AGain*acc_sigma2;
      MatrixScalarMul(tgain,Qkk);
      MatrixMul(Hk,xkk1,temp11);
      //MatrixPrint(xkk1);
      MatrixAdd(zk,temp11,yk, -1); 
      MatrixMul(Hk,Pkk1,temp12);
      MatrixMul(temp12,Hk_t,temp11);
      MatrixAdd(temp11,Rk,Sk,1);
      MatrixInvert(Sk,Sk_inv);
      MatrixMul(Pkk1,Hk_t,temp21);
      MatrixMul(temp21,Sk_inv,Kk);
      MatrixMul(Kk,yk,temp21);
      MatrixAdd(temp21,xkk1,xkk,1);
      
      MatrixMul(Kk,Hk,temp22);
      MatrixAdd(Eye,temp22,temp22a,-1);
      MatrixMul(temp22a,Pkk1,Pkk);
      
      //predict cycle
          //make copies of the last iteration;
      ArrayCopy(Pk1k1,Pkk);
      ArrayCopy(xk1k1,xkk);
      //g = MathSqrt(MathPow(xkk[0][0]-get_avg(i),2))/Pkk[0][0]*Point;
      //g = xkk[1][0];
      //predict the state
      //Print("Doing Predict");
      MatrixMul(Fk,xk1k1,xkk1);
      //MatrixPrint(xkk1);
      MatrixMul(Fk,Pk1k1,temp22);
      MatrixMul(temp22,Fk_t,Pk1k1);
      MatrixAdd(Pk1k1,Qkk,Pkk1,1);
      //Print("Predict Ended");
      KalmanBuffer[i-1] = (xkk1[0][0]*Point);
 
      if(i==1){
         if(PrintCurrentState){
            Print("Suppression (dB):", -20*MathLog(MathSqrt(tgain)/Point)/MathLog(10.),"dB; Value:",xkk1[0][0]*Point
            ,"; Speed (pips/bar):",xkk1[1][0],"; dp/p:"
            , xkk1[1][0]/xkk1[0][0]*100.*240.*(24.*60.)/Period());
         }
         if(PredictBars1 > 0){
            MatrixMul(Phi1,xkk,xpt);
            xp[0][0] = xkk[0][0];
         
            MatrixTranspose(Phi1,temp21);
            MatrixMul(Pkk,temp21,temp21a);
            MatrixMul(Phi1,temp21a,temp11a);
            MatrixInvert(temp11a,temp11);
            MatrixAdd(xp,xpt,xdel,-1);
            MatrixTranspose(xdel,xdel_t);
            MatrixMul(temp11,xdel,temp11a);
            MatrixMul(xdel_t,temp11a,temp11);
            g1=MathSqrt(temp11[0][0]);
            Print(Phi1[0][1], " bar bounds ",(xpt[0][0]-g1)*Point
               ,",",(xpt[0][0]+g1)*Point);
            Print(Phi1[0][1], " bar center, r(t) (in pips) ",(xpt[0][0])*Point
               ,",",g1);
         }
         if(PredictBars2 > 0){
            MatrixMul(Phi2,xkk,xpt);
            xp[0][0] = xkk[0][0];
            MatrixTranspose(Phi2,temp21);
            MatrixMul(Pkk,temp21,temp21a);
            MatrixMul(Phi2,temp21a,temp11a);
            MatrixInvert(temp11a,temp11);
            MatrixAdd(xp,xpt,xdel,-1);
            MatrixTranspose(xdel,xdel_t);
            MatrixMul(temp11,xdel,temp11a);
            MatrixMul(xdel_t,temp11a,temp11);
            g1=MathSqrt(temp11[0][0]);
            Print(Phi2[0][1], " bar bounds ",(xpt[0][0]-g1)*Point
               ,",",(xpt[0][0]+g1)*Point);
            Print(Phi2[0][1], " bar center, r(t) (in pips) ",(xpt[0][0])*Point
               ,",",g1);
         }         
      }else g= -1;
      
   
   }
   //Print("24 hours:", xpt[0][0]*Point," ,Gate ",g);
   //Print ("Gate 2 ", MathSqrt(xdel[0][0]*xdel[0][0]*g)*Point," g1 ",g1);

}
 
 double get_avg(int k){
   return(MathPow((High[k]*Low[k]*Close[k]*Close[k]),1/4.)/Point);
}        
//+++++++-----------------------------------------------------+