どんなもの？
解りやすくまとめてあるサイトが数多くありますので、参照してください。手順は、以下の様になります。
- モデルの決定
- 粒子の初期化（散布）
- 予測
- 尤度の計算、重み付
- リサンプリング
- 重み付平均
- 3に戻る
作成
今回作成したプログラムの基本設計は、以下のとおりです。また、直前の値からほとんど動かない予測値を想定しています。【モデル】
| システムモデル | x[t]=x[-t1]+v （※vは標準正規分布のノイズ） |
| 観測モデル | y[t]=x[t]+w （※w=0と設定) |
| 粒子設計 | Particles[粒子数][構造] （構造[0：位置情報,1：重み情報]) |
| 尤度の計算 | 標準正規分布を使用 |
| 正規分布の偏差 | 直近データより取得 |
| リターン値 | 直近の終値を基準にした予測値（黄線）と答えとなる終値（赤線） |
【パラメータ】
| int Particles_vol =1000; | 粒子数量を設定するパラメータ |
| int Range = 100; | 試行範囲（粒子数量×試行回数で計算量が増えるので注意） |
| int sigmaperiod =100; | 標準正規分布のパラメター（偏差）を取得するための計算範囲 |
| double Set = 3000.0; | 粒子の初期設定の一様分布の範囲（pipsで記載） -Set<= Particles[][0] <=Set |
【コード】
//+------------------------------------------------------------------+ //| SMCFilter.mq4 | //| Copyright 2012, MetaQuotes Software Corp. | //| http://www.metaquotes.net | //+------------------------------------------------------------------+ #property copyright "Copyright 2012, MetaQuotes Software Corp." #property link "http://www.metaquotes.net" #property indicator_separate_window //#property indicator_chart_window #define PI 3.1415926 #property indicator_buffers 2 #property indicator_color1 Yellow #property indicator_width1 1 #property indicator_color2 Red #property indicator_width2 1 extern int Particles_vol = 1000;//粒子の数 extern int Range = 100;//試行範囲 extern int sigmaperiod = 100;//偏差取得期間 extern double Set = 3000.0;//初期化範囲（pips) double Linebuf[];double ans[];double Particles[][2];//+------------------------------------------------------------------+ //| Custom indicator initialization function | //+------------------------------------------------------------------+ int init() { //---- indicators SetIndexStyle(0, DRAW_LINE); SetIndexBuffer(0,Linebuf); SetIndexEmptyValue( 0,0); SetIndexStyle(1, DRAW_LINE); SetIndexBuffer(1,ans); SetIndexEmptyValue( 1,0); //乱数の初期化 MathSrand(TimeLocal()); ArrayResize(Particles,Particles_vol); //---- return(0); } //+------------------------------------------------------------------+ //|粒子の初期化 | //+------------------------------------------------------------------+ void setparticles(double & pobj[][],int range,double set){ for(int i=0;i<Particles_vol;i++){ Particles[i][0] = data(range)+(((set*2)/Particles_vol)*i-set)*Point; Particles[i][1] = 0; } } //+------------------------------------------------------------------+ //|乱数の成形 0<=random<1 | //+------------------------------------------------------------------+ double random(){ double rd = MathRand(); rd/= 32768; return(rd); } //+------------------------------------------------------------------+ //|正規分布を計算 mu:平均 sigma:偏差 | //+------------------------------------------------------------------+ double normal(double x,double sigma){ double mu = 0; return(MathExp(-MathPow((x-mu),2)/(2*sigma*sigma))/(MathSqrt(2*PI)*sigma)); } //+------------------------------------------------------------------+ //|正規分布となる乱数の成形 mu:平均 sigma:偏差 | //+------------------------------------------------------------------+ double normalrand(double sigma){ double mu = 0; double t,u,r1,r2; sigma = MathSqrt(sigma); t=MathSqrt(-2.0*MathLog(1-random())); u=2*PI*random(); r1=t*MathCos(u); return(r1*sigma+mu); } //+------------------------------------------------------------------+ //|リサンプリング q:粒子数 | //+------------------------------------------------------------------+ void resample(int q,double &pobj[][] ){ double wt[]; double rd; int i,j; //粒子のコピー double copyp[][2]; ArrayResize(copyp,q); ArrayCopy(copyp,pobj,0,0,q*2); //累計重み ArrayResize(wt,q); wt[0] = pobj[0][2]; for(i=1;i<q;i++) wt[i]=wt[i-1]+pobj[i][2]; //リサンプリング(ルーレット法） for(i=0;i<q;i++){ rd = random()/Particles_vol; for(j=0;j<q;j++){ if(rd > wt[j]){ continue; }else{ pobj[i][0] = copyp[j][0]; pobj[i][1] = 0; break; } } } } //+------------------------------------------------------------------+ //|尤度計算および重み付け | //+------------------------------------------------------------------+ void wlikelihood(double price,int p,double sigma,double &pobj[][]){ double alfa[]; double wsum; double g = 1; //Gの偏微分の値 int i; ArrayResize(alfa,p); for(i=0;i<p;i++){ alfa[i] = normal((price-pobj[i][0]),sigma)*g;//(尤度） wsum += alfa[i]; } if(wsum<1)wsum=1; //正規化 for(i=0;i<p;i++)pobj[i][1] = alfa[i]/wsum; } //+------------------------------------------------------------------+ //|対象データ | //+------------------------------------------------------------------+ double data(int index){ return(Close[index]); } //+------------------------------------------------------------------+ //|偏差の取得 | //+------------------------------------------------------------------+ double setsigma(int range,int index){ int i; double sum,av; for(i=0;i<range;i++)av += data(index+i); av /= range; for(i=0;i<range;i++)sum += MathPow((data(index+i)-av),2); return(MathSqrt(sum/range)); } //+------------------------------------------------------------------+ //| Custom indicator iteration function | //+------------------------------------------------------------------+ int start(){ double v;//乱数 double sum,sigma,limit; int i,j; static datetime now = 0; int counted_bars=IndicatorCounted(); if(counted_bars<1)i=Range; else i=0; if(now==0)setparticles(Particles,Range,Set); //limit =Set*Point; while(i>=0){ //各bar一回のみ実行 if(now!=Time[i]){ now=Time[i]; sigma=setsigma(sigmaperiod,i); for(j=0;j<Particles_vol;j++){ Particles[j][0] = Particles[j][0] + normalrand(sigma); } //尤度推定 wlikelihood(data(i+1),Particles_vol,sigma,Particles); //リサンプリング resample(Particles_vol,Particles); //加重付平均 sum=0.0; for(j=0;j<Particles_vol;j++)sum += (Particles[j][0])*Particles[j][1]; Linebuf[i] = sum; } ans[i] = data(i);//答え i--; } return(0); } //+------------------------------------------------------------------+
テスト
稼働させると以下の様になります。
今回のテストでは、直前の終値から次の終値は、ほとんど変化しないという条件で作成したため、このような結果となりました。
0 件のコメント :
コメントを投稿