Contoh dan Kode Runge-Kutta

Melanjutkan postingan terakhir saya, saya pikir saya akan, sebagai langkah pertama, membuat kode “sederhana” Runge Kutta fungsi dan menunjukkan cara menangani fakta bahwa tidak ada “rumus matematika ajaib” untuk menghitung kemiringan yang merupakan bagian integral dari Runge-Kutta.

Pendekatan saya adalah menyesuaikan fungsi kuadrat ke n_bar terakhir harga dan ambil kemiringan ini melalui saya Konvolusi filter Savitzky-Golay kode, dan dengan demikian nilai Runge-Kutta k1 dapat diperoleh dengan mudah. ekstrapolasi melampaui titik k1 ke titik k2 dan k3, pada dasarnya, mencoba menyesuaikan ke titik yang memiliki keunggulan setengah bar atas harga terakhir yang tersedia. Untuk mengakomodasi ini, saya menggunakan n_bar – 1 poin terakhir dari simple moving average 2 bar plus “posisi” titik k2 dan k3 untuk menghitung kemiringan pada k2 dan k3. Rata-rata pergerakan sederhana 2 bar digunakan karena ini memiliki jeda setengah bar dan secara efektif merupakan interpolasi dari harga yang diketahui pada “interval setengah bar,” dan oleh karena itu titik k2 dan k3 adalah satu langkah h di depan interval setengah bar terakhir. Titik k4 lagi-lagi dihitung langsung dari harga plus proyeksi interval setengah bar dari titik k3. Jika semua ini tampak membingungkan, mudah-mudahan Oktaf Kode di bawah akan memperjelas semuanya.

clear all

%  create the raw price series
period = input( 'Period? ' ) ;
sideways = zeros( 1 , 2*period ) ;
uwr = ( 1 : 1 : 2*period ) .* 1.333 / period ; uwr_end = uwr(end) ;
unr = ( 1 : 1 : 2*period ) .* 4 / period ; unr_end = unr(end) ;
dwr = ( 1 : 1 : 2*period ) .* -1.333 / period ; dwr_end = dwr(end) ;
dnr = ( 1 : 1 : 2*period ) .* -4 / period ; dnr_end = dnr(end) ;

trends = [ sideways , uwr , unr.+uwr_end , sideways.+uwr_end.+unr_end , dnr.+uwr_end.+unr_end , dwr.+uwr_end.+unr_end.+dnr_end , sideways ] .+ 2 ;
noise = randn( 1 , length(trends) ) .* 0.0 ;
price = sinewave( length( trends ) , period ) .+ trends .+ noise ;
ma_2 = sma( price , 2 ) ;
 
%  regress over 'n_bar' bars
n_bar = 9 ;
%  and a 'p' order fit
p = 2 ;

%  get the relevant coefficients
slope_coeffs = generalised_sgolay_filter_coeffs( n_bar , p , 1 ) ; 

%  container for 1 bar ahead projection
projection_1_bar = zeros( 1 , length( price ) ) ;

for ii = n_bar : length( price )

%  calculate k1 value i.e. values at price(ii), the most recent price
k1 = price( ii-(n_bar-1) : ii ) * slope_coeffs( : , end ) ;
projection_of_point_k2 = price(ii) + k1 / 2 ;

%  calculate k2 value
k2 = [ ma_2( ii-(n_bar-2) : ii ) ; projection_of_point_k2 ]' * slope_coeffs( : , end ) ;
projection_of_point_k3 = price(ii) + k2 / 2 ;

%  calculate k3 value
k3 = [ ma_2( ii-(n_bar-2) : ii ) ; projection_of_point_k3 ]' * slope_coeffs( : , end ) ;
projection_of_point_k4 = price(ii) + k3 / 2 ;

%  calculate k4 value
k4 = [ price( ii-(n_bar-2) : ii ) , projection_of_point_k4 ] * slope_coeffs( : , end ) ;

%  the runge-kutta weighted moving average
projection_1_bar(ii) = price(ii) + ( k1 + 2 * ( k2 + k3 ) + k4 ) / 6 ;

end

%  shift for plotting
projection_1_bar = shift( projection_1_bar , 1 ) ;
projection_1_bar( : , 1:n_bar ) = price( : , 1:n_bar ) ;

plot( price , 'c' , projection_1_bar , 'r' ) ;

Kode ini menghasilkan plot seperti ini, tanpa noise,

dan ini dengan noise (baris 12 kode).

Garis cyan adalah harga dasar dan garis merah adalah proyeksi Runge-Kutta 1 bar ke depan. Seperti yang dapat dilihat, ketika harga bergerak dalam garis yang agak lurus, proyeksi tersebut cukup akurat, namun, pada saat pembalikan ada beberapa overshoot, yang memang diharapkan. Saya tidak terlalu khawatir tentang overshoot ini karena tujuan saya hanya untuk mendapatkan berbagai nilai k sebagai fitur, tetapi overshoot ini memang memiliki beberapa implikasi yang akan saya bahas di posting mendatang.