移動平均と標準偏差の算出 ~株価チャート:シグマバンド(ボリンジャーバンド)の作成~

コメントをどうぞ


「Igor:waveの読み込み3」で扱ったデータを元に,株価チャートの一つであるシグマバンド(ボリンジャーバンド)を作成する。

シグマバンドとは,株価の移動平均に対して,移動平均算出区間内での標準偏差σの線を引いて,株価の異常値を見つけるチャートである。

株価の変動について正規分布を仮定すると,平均値±σの範囲には68%±2σの範囲には95.5%±3σの範囲には99.7%の確率で株価は収まることになる。

 

移動平均を表すwave,+2σおよび-2σを表すwaveを作成する。
「Igor:waveの読み込み3」で読み込んだデータ数は250であったので,「Make」コマンドで,

Make /O /N=250 ‘wave_mean’ ‘wave_mean_2stdp’ ‘wave_mean_2stdm’;
AppendToTable wave_mean wave_mean_2stdm wave_mean_2stdp

とする。

移動平均算出のために,以下のマクロ「calc_mean(input_w, output_w)」を作成した。
input_wの移動平均値がoutput_wに代入される。
input_wとして株価の終値であるclosing_priceを代入し,wave_meanにその移動平均値が代入されるように使用する。

平均を算出する関数「 mean(waveName, [x1, x2]) 」を利用した。
waveNameのx1~x2の区間の平均値を算出してくれる関数である。

waveのデータ数が250であり,移動平均算出区間は20データ(定数spanとした)であるため,for文の区間は230とすべきである。例えばi=240の場合,mean(input_w, 240, 260)となり,waveデータが存在しないところを参照してしまう。
ただし,Igorはこのような場合,自動的にmean(input_w, 240, 249)として出力してくれるようである。

標準偏差を算出するために,以下のマクロ「calc_std(input_w, output_w1, output_w2, mean_w)」を作成した。
input_wの標準偏差σの2倍を移動平均値にプラスしたoutput_w1,マイナスしたoutput_w2を算出する。

標準偏差σを算出にあたり,分散σ2を算出する「 Variance(waveName, [x1,x2]) 」を利用した。
waveNameのx1~x2の区間の分散を算出してくれる関数である。

分散の平方根をとることで(sqrt関数),標準偏差を算出した。

移動平均と標準偏差

移動平均と標準偏差

 

以下のように,マクロを実効

calc_mean(closing_price,wave_mean)
calc_std(closing_price,wave_mean_2stdp,wave_mean_2stdm,wave_mean)

メニューバーの[Window]->[New Graph]よりグラフを作成。
XWave:[Day],YWave:closing_price, mean_wave, mean_2stdp, mean_2stdm
を選択。

見栄えは,軸をダブルクリックしていろいろと変更,ラインをダブルクリックしていろいろと変更。
メニューバーの[Graph]->[Add Anotation]から凡例を設定。

これらどのように設定したかは,以下のコマンドラインを参照のこと。

シグマバンドのできあがり

シグマバンドのできあがり

広告

Igor : waveの読み込み3

コメントをどうぞ


Igorでwaveの読み込み方法について記す。前2回も参照のこと。

今回は,ヘッダ部をスキップして読み込むことと,日付データの読み込みを行う。

例として,株価データの読み込みを行う。以下のサイトから個別銘柄の株価データをダウンロードできる。
http://k-db.com/stocks/

直近250日分の日足をcsvファイル形式でダウンロード。

日足データのcsvファイル

日足データ

このcsvファイルは,最初の2行がヘッダ情報になっており,3行目からが日足データの本体。
1行目がタイトル行であり,2行目が各列の列名。

読み込みには,「LoadWave」を用い,

LoadWave /J /E=2 /L={1,2,0,0,0} /R={English,2,2,2,1,”Year-Month-DayOfMonth”,40}

とする。

/Jと/Eフラグは前回も用いたが,ASCIIデータであることと既存のテーブルに追加することを明示している。

/Lでヘッダをスキップすることを明示している。
/L={nameLine, firstLine, numlines, firstColumn, numColumns}という形式で指定する。
日本語でざっくり書くと,
/L={名前行は何行目か, データを何行目から読み込むか, データの行数, データを何列目から読み込むか, データの列数}
となる。

ここで,「0」とするとAutoとなる。numlines(行数)など,指定せずファイルにあるだけ読み込みたいときは0にすればよい。

/L={1,2,0,0,0}としているので,名前行は1行目(0からスタート)。ただしこのファイルは日本語データであり,うまく読み込めないので,0としてもよい。

 

/Rは日付形式を指定する。このcsvデータ(1列目)はIgorでの標準の日付形式でないため,指定してあげる必要がある。
/R={languageName, yearFormat, monthFormat, dayOfMonthFormat, dayOfWeekFormat, layoutStr, pivotYear}という形式で指定する。
今回は,/R={English,2,2,2,1,”Year-Month-DayOfMonth”,40}と指定した。

yearFormat: 1:数字2つで表記,2:数字4つで表記  14 or 2014
monthFormat: 1:数字(0省略),2:数字(0省略なし),3:アルファベット(省略表記,Janなど),4:アルファベット(省略なし,Januaryなど)
dayOfMonthFormat: 1:数字(0省略),2:数字(0省略なし)
dayOfWeekFormat: 1:アルファベット(省略表記,Monなど), 2:アルファベット(省略なし,Mondayなど)
layoutStr: レイアウト例を明示(図参照)
pivotYear: この数字(yyとする)以下ならば,20yy年。以上ならば19yy年。yyは4位上

Rフラグの説明

Rフラグの使い方

Rフラグの例

Rフラグの例

CSV読み込み後のIgor画面は下図。wave0~wave6まで自動で名前付けされた。

CSV読み込み後

CSV読み込み後

読み込み時のダイアログでもwave名は編集できるが,以下のようにRenameコマンドで名称変更してもよい。

Rename wave0,Day; Rename wave1,Opening_price; Rename wave2, high_p; Rename wave3, low_p;Rename wave4,closing_price; Rename wave5,trading_volume; Rename wave6 trading_value;

Igorでバイナリファイルの読み込み

コメントをどうぞ


Igorでバイナリファルを読み込む。

ファイル名はabc.binでファイルの中身はdouble型が50個格納されている。double型は8バイトなので合計400バイトとなる。

メニューで
Data/Load Waves/Load General Binary File
をクリック。

Load General Binary Data

型:Double float
Byte order:Low byte first (intel系の場合)
File:読み込むファイルを指定(今回はabc.bin)

実行後(Do it)
AppendToTable コマンドでテーブルに追加。以下のように読み込めた。

Output of Load Binary Data

補足1:

データが行列のように1次元でない場合,その列毎にwaveを分けたかったりする。
その場合には,以下のようにプロシージャを作成すれば良い。
GBLoadWaveコマンドを,読み込むオフセットを設定し,for文で繰り返し実行している。

Function test()
String file_name=”C:test:abc.bin”
Variable i, offset_i
offset_i = 0
for (i=0;i<5;i+=1)
GBLoadWave /B /T={4,4} /N=$(“wave”+num2str(i)) /S=(offset_i) /W=1 /U=10 file_name
offset_i += 50
endfor
end

補足2:

今回は,C言語でバイナリファイルを作成した。ファイルの中身は,double型で10行5列に乱数を代入している。

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#pragma warning(disable:4996)

int main(void)
{
int i,j;
double abc[10][5];

FILE *file;

for(i=0;i<10;i++){
for(j=0;j<5;j++){
abc[i][j]=rand();
printf(“%lf\n”,abc[i][j]);
}
}

file = fopen(“abc.bin”,”wb”);
fwrite(abc,sizeof(abc),1,file);
fclose(file);

return 0;
}

Igorで行列計算(ジョーンズ行列で偏光状態を計算)

コメントをどうぞ


前稿で,Igorで行列計算をする方法について述べた。
本稿では,その例題として,ジョーンズ行列を用いた,波長板による偏光の計算を行う。
ジョーンズ行列は,光学系の偏光状態を計算するのに用いられる。

直線偏光を表すジョーンズベクトルは,

左円偏光を表すジョーンズベクトルは,

右円偏光を表すジョーンズベクトルは,

のように表される。

偏光状態を変化させる光学素子として,波長板がある。波長板は,速軸(Fast axis)と遅延軸(Slow axis)を有する素子である。波長板は,それを透過した光の遅延軸成分が速軸に対して,位相がΓ遅れる性質を持っている。
この波長板のジョーンズマトリックスTは,

で表される。

Γ=π/2のとき,1/4波長板と呼ばれる。
直線偏光J1 (1;1)を左円偏光J2 (1;-j)に転換し,右円偏光J3 (1;j)を直線偏光J1 (1;1)に転換する。

Γ=πのとき,1/2波長板と呼ばれる。
直線偏光J1 (1;1)を直線偏光J4 (1;-1)に転換し,右円偏光J3 (1;j)を左円偏光J2 (1;-j)に転換する。

1/4波長板と1/2波長板の例

1/4波長板と1/2波長板の例

1/4波長板と1/2波長板の計算をIgorで行ってみた。

1/4波長板のジョーンズ行列をT1,1/2波長板のジョーンズ行列をT2とすると,それぞれ,

となる。

1/4波長板に直線偏光J1を入射すると,左円偏光J2になる。

1/2波長板に右円偏光J3を入射すると,左円偏光J2になる。

以上をIgorで計算すると,以下のようになる。

//ジョーンズ行列

Make /N=(2,2) /C T1,T2
Make /N=2 /C J1,J3
variable theta1=pi/2 //1/4波長板
variable theta2=pi //1/2波長板
J1={1,1}
J3={1,cmplx(0,1)}
T1[][0]={1,0}
T1[][1]={0,cmplx(cos(theta1),-sin(theta1))}
T2[][0]={1,0}
T2[][1]={0,cmplx(cos(theta2),-sin(theta2))}
MatrixOp /O J2=T1 x J1
MatrixOp /O J2_=T2 x J3
appendtotable T1,T2,J2,J2_

実行例は下図。

Igorでジョーンズ行列計算

Igorでジョーンズ行列計算

Igorで行列計算

コメントをどうぞ


Igorで行列計算

例として以下の行列T1とT2の掛け算をIgorで計算する。

結果は,

となる。

以上の計算をigorで行った。

以下はそのプログラム例

Make /N=(2,2) T1,T2,T3 //2行2列の行列T1,T2,T3を定義

T1[][0] = {1,2} //0列目に値代入(行列番号は0から始まる)
T1[][1] = {3,4} //1列目に値代入
T2[][0] = {5,6}
T2[][1] = {7,8}

MatrixMultiply T2,T1 //行列の掛け算
T3 = M_product //結果はM_productとして保存される。

AppendToTable T1,T2,T3

matrixmultiplyの計算例

行列計算結果

ただし行列の掛け算部は,以下の演算子”MatrixOp“を用いたほうが簡単。
行列掛け算部の2行は次の1行と等価である。

MatrixOp /O T3 = T2 x T1 //行列の掛け算はx(エックス)を用いる。

MatrixOpを用いる場合には,T3を前に宣言しておく必要なし。
しかし,ここでは逆に,宣言してしまっているから”/O”フラグで上書き許可してる。

Igorで複数軸のグラフ作成

コメントをどうぞ


Igorで複数軸のグラフを作成する方法について記す。

2つの系列でそれぞれ左軸,右軸というグラフはよく作るが,今回は,左軸だけで分離して作成する方法を記す。

 

例として,”Time_s”,”V1″,”V2″という3つのデータ系列を用意し,それをグラフ化する。

“Windows”タブの”New Graph”を選択すると,下図の画面が現れる。”Fewer Choices”ボタンを押すと,簡易設定画面になるが,今回は下図の画面のように詳細設定できる画面ですすめる。

グラフ作成画面

XWaveとしてTime_sを選択し,YwaveとしてV1を選択し,”Add”ボタンを押す。

第2系列を追加

 

同様に,XWaveとしてTime_sを選択し,YwaveとしてV2を選択し,”Add”ボタンを押す。
V2の系列のところで,YAxisの”left”となっているところをクリックし,”New”をクリックする。

そして,Nameの欄に適当に名前をつける。ここでは”L2″とした。

新規軸を追加

 

続いて,”Do It”ボタンを押すと,下図のようなグラフが作成される。

グラフ描画

このままでは,左軸の2つが重なりあってしまい,見栄えが良くないため,修正をする。

左軸をダブルクリック(右クリックして”Axis Properties”でもよい)する。

Axisでleftを選択し,Axisタブの中の”Draw between”のところを0~70%とする。

左軸の修正

 

続いて,Axisで”L2″を選択し,Axisタブの中の”Draw between”のところを80~100%とする。

L2というのは先ほどつけた新規左軸の名前である。

そして,”Distance from Margin”のところを0とする。こうすることで,既存のleft軸と同じ直線上に新規左軸が移動する。

新規左軸の修正

複数軸のグラフはこれで完成である。

あとは,適当にグラフの見栄えを修正すると,下図のようなグラフとなる。

グラフの完成

Igorで複素数計算 ローパスフィルタの特性計算

コメントをどうぞ


Igorで複素数の計算をする方法について書く。

Igorでのコマンドは以下のリンクに詳しい。

http://www.hulinks.co.jp/support/igor/reference/functions/index_complex.html

 

例題として,図のようなローパスフィルタの特性について計算した。

R=100Ω,C=1μFである。

RとCによるローパスフィルタの回路図

RとCによるローパスフィルタの回路図

入力と出力の電圧比は,次式となる。

入出力の電圧比

入出力の電圧比

 

この絶対値が電圧利得の周波数特性になり,極座標表示における,虚部が位相のずれの周波数特性を表す。

Igorでのプログラム例を示す。

用いた関数は,cmplxr2polarである。

cmplx(実部,虚部)

という書式で,複素数を定義する。

wave1=r2polar(wave2)

という書式で,wave2を極座標変換し,wave1を作成する。wave1の実部がwave2の大きさ,wave1の虚部がwave2の位相である。

また,複素数の大きさについては,

wave3=magsqr(wave4)

というようにすればwave4(複素数)の大きさの二乗がwave3(実数)に代入される。

複素数のwave作成時には,Makeコマンドに/cというフラグを付ければよい。

 

プロシージャ画面で,関数(Lowpass)を作成。

ローパスフィルタ特性を計算する関数の定義

ローパスフィルタ特性を計算する関数の定義

実行すると,

関数実行後の特性グラフ化

関数実行後の特性グラフ化

となる。

電圧利得と位相ずれの周波数特性が計算できた。

コマンドを以下に示す。

Function Lowpass()
    Variable R=100
    Variable C=1e-6
    Make /o /N=100 wave_freq = 0.05*x    //0から0.05刻みでwaveの初期値を作成
    wave_freq = 10^(wave_freq)            //logスケールで周波数を定義
    Make /o /N=100 wave_omega = 2*pi*wave_freq
    Make /o /N=100 /c wave_BoardDiag    //入出力の電圧比(複素数)
    Make /o /N=100 mag_BoardDiag        //電圧利得の周波数特性
    Make /o /N=100 phase_BoardDiag        //位相ずれの周波数特性
    wave_BoardDiag = cmplx(1,wave_omega*C*R)    //1+jωC
    wave_BoardDiag = 1/wave_BoardDiag            //1/(1+jωC)
    mag_BoardDiag = real(r2polar(wave_BoardDiag))    //極座標の実部が電圧利得を表す。
    //mag_BoardDiag = sqrt(magsqr(wave_BoardDiag))としても電圧利得が計算できる
    phase_BoardDiag = 180/pi*imag(r2polar(wave_BoardDiag))    //極座標の虚部が位相ずれを表す。
    AppendToTable wave_freq, wave_BoardDiag, mag_BoardDiag,phase_BoardDiag    //各waveの値をテーブルに表示
    Display mag_BoardDiag vs wave_freq                    //グラフ作成
    ModifyGraph grid=1,log=1,tick=2,mirror=1   
    ModifyGraph width={Aspect,1.2}                       
    Label left “V\\Bout\\M/V\\Bin”
    Label bottom “Frequency [Hz]”
    Display phase_BoardDiag vs wave_freq     //グラフ作成
    ModifyGraph grid=1,log(bottom)=1,tick=2,mirror=1
    ModifyGraph width={Aspect,1.2}
    Label left “Phase [deg]”
    Label bottom “Frequency [Hz]”
End

 

Older Entries