このページでは、財政学の研究者向けにフォートラン言語の解説と私が作成したフォートラン・プログラムの公開をおこないます。最近では、パッケージソフトでもシミュレーション分析は十分可能ですが、プログラム言語ならでは自由度の高さはやはり魅力です。以下の解説を見るときは、手元にfortranの解説書をおいてくださいね。すべての命令を細かく説明するのは面倒ですから・・・
パソコンではwindows95 ないし ＮＴ 上で動く、フォートランとして、かってはマイクロソフトが英語版でFortran
Powerstation Ver 4.0を発売していました。これは、Fortran90に対応した３２ビットアプリケーションで専用エディターも付属した優れものでした。Fortran90はそれまでの規格だったFortran77に比べると、行番号なしにループ文が使えるなどよりBasicライクになり、簡単になりました。しかし、残念ながらマイクロソフトはフォートランから撤退し、その後digital
equipmentにサポートが引き継がれることになりました。
私が現在使用しているのは、このdigital fortranです。マイクロソフト版とほぼ使い勝手は同じでかつfortran95の規格による命令の一部が使えるようです（私はfortran90の命令しか使いませんが・・・）。
さらにコンパックによるDECの買収？により、いまから買うならコンパック版になります。
フォートランでは、合計を出すのにループ文を使います。たとえば３＋１０＋５を求めるなら
REAL X（1:3） ! X（）という配列変数を実数型で宣言します。
REAL total
INTEGER I ！整数型変数Iを宣言
X(1)=3 ! 変数にデータを代入します
X(2)=10 ! 変数にデータを代入します
X(3)=5 ! 変数にデータを代入します
Total=0 ! 合計の変数としてのTOTALをゼロクリアーします
Do I=1,3
total=total+X(I)
write(*,*) I,total,X(I)
end do
stop
end
で計算できます。ただし、配列変数へのデータセットは、この方法だとデータ数が多くなると大変です。もう少しスマートにやるなら
REAL X（1:3） ! X（）という配列変数を実数型で宣言します。
REAL total
INTEGER I ！整数型変数Iを宣言
DATA X /3,10,5/
Total=0 ! 合計の変数としてのTOTALをゼロクリアーします
Do I=1,3
total=total+X(I)
write(*,*) I,total,X(I)
END DO
STOP
END
という感じですね・・・
シミュレーション・プログラムを書くのにループ文以外に書かせないのがif文です。もし××ならば、○○せよ。 という命令です。 ××のところには、basicだと＝とか＞を使えるのですがフォートランの場合、EQとかＧＴを使います。要するに英語なんですよね・・・。もし日本人がフォートランを設計したならば、等しい、とかより大きいとか なったんでしょうか・・・。
たとえば、xの値がゼロになったら、ストップさせるプログラムを書くなら
ＩＦ (X.EQ.0) STOP
でよかったはずです・・・
しかし、実際のIF文の使い方としては、もし××ならば、○○と△△せよ、というふうに複数の計算が要求される 場合の方が多いと思います。そんなときには IF と THEN とEND IF が役に立ちます。 ちょっと実例を参考にしてみましょう。以下のプログラムは課税所得が与えられた場合に 平成10年の税率表にもとづいて所得税額を計算するプログラムです。 かなりややこしいがループ文とIF文だけで構成されてます。 電卓で計算するなら課税所得が1000万円なら330×0.1+(900-330)×0.2+(1000-900)×0.3=177万円 と求めるのですが、フォートランだと以下のようになります。
REAL TAXBASE,TAX,B(1:6),T(1:5),TOTAL,MT
INTEGER I
TAXBASE=1000 ！課税所得に1000万円を代入
DATA B /0,330,900,1800,3000,1000000000/ !税率表の刻みです
DATA T /0.1,0.2,0.3,0.4,0.5/ !限界税率です
TOTAL=0
DO I=1,5
IF ((B(I).LT.TAXBASE).AND.(TAXBASE.LE.B(I+1)))
THEN
TAX=TOTAL+(TAXBASE-B(I))*T(I)
MT=T(I)
END IF
TOTAL=TOTAL+(B(I+1)-B(I))*T(I)
END DO
WRITE(*,*) TAX,MT
STOP
END
IF ((B(I).LT.TAXBASE).AND.(TAXBASE.LE.B(I+1)))
THEN
TAX=TOTAL+(TAXBASE-B(I))*T(I)
MT=T(I)
END IF
がIF文になってますよね。たとえばIが１のときこの命令は
もし、 ０＜課税ペース かつ 課税ベース＜＝330 のときに
TAX=TOTAL+(TAXBASE-B(I))*T(I)
MT=T(I)
を実行しますが、課税ベースは1000でこの条件にあてはまらないのでTHENとEND IFに囲まれた２行の命令は無視されます。同じくIが２のときも条件が満たされないので間にはさまれた命令は実行されません。課税ベースが1000のなので、Iが３のときだけ条件が満たされます。つまり
もし ９００＜課税ベース かつ 課税ベース＜＝１０００のとき
税額=１４７＋（1000-900）*0.3
ＭＴ＝0.3
が実行されます。ここでＩが３になったときにはTOTALという変数に１４７がセットされているはずですがなぜだかわかりますか。合計を出してみようを見ながら自分で考えてみてください。
fortranと昔のN88BASIC（知っているあなたは30歳以上のはず・・・）の最大の違いは、 サブルーチンで引数を使えるかどうかでしょう。サブルーチンを使えば、プログラムをブロックごとの部品化して、汎用性をもたせることができます。
論より証拠ということで、私が作成した所得税額を計算するプログラムで解説しましょう。
!平成９年税制
INTEGER E,S,B,K,C,I
PARAMETER(E=5,S=2,B=5,C=3,K=2)
real WB(6),WR(5),PB(3),PR(2),A(6),T(5)
real WAGE,INCOME,WAGEDE,PREMI,LKY,MEN,TAX,TBASE
real KISO,HAIGU,FUYOU,FMEN
DATA (WB(I),I=1,6) /0,180,360,660,1000,100000000/
DATA (WR(I),I=1,5) /0.4,0.3,0.2,0.1,0.05/
DATA (A(I),I=1,6) /0,330,900,1800,3000,1000000000/
DATA (T(I),I=1,5) /0.1,0.2,0.3,0.4,0.5/
DATA (PR(I),I=1,2) /0.07,0.02/
DATA (PB(I),I=1,3) /0,500,10000000/
wage=1000
men=4
TAX=0
WAGEDE=0
LKY=65
CALL RATE(E,WAGE,WAGEDE,WB,WR)
IF (WAGEDE.LT.LKY) WAGEDE=LKY
INCOME=WAGE-WAGEDE
IF (INCOME.LT.0) INCOME=0
KISO=38
HAIGU=38
FMEN=MEN-2
IF (FMEN.LT.0) FMEN=0
FUYOU = 38*FMEN
IF (MEN.LT.2) HAIGU = 0
IF (INCOME.LE.1000) HAIGU=HAIGU*2
CALL RATE(S,WAGE,PREMI,PB,PR)
IF (WAGE.GT.1000) PREMI=45
TBASE=INCOME-(KISO + HAIGU + PREMI + FUYOU)
CALL RATE(B,TBASE,TAX,A,T)
write(*,*) TAX
END
SUBROUTINE RATE (E, TBASE, ANS, B, R)
INTEGER E,I
real TT,B(1:E+1),TBASE,ANS,R(1:E)
TT = 0
DO I = 1,E
IF ((B(I).LT.TBASE).AND.(TBASE.LE.B(I + 1)))
THEN
ANS = TT + (TBASE - B(I)) * R(I)
END IF
TT = TT + (B(I + 1) - B(I)) * R(I)
END DO
RETURN
END
今回は長くて全部解説するのは面倒なので、ポイントだけね・・・
サブルーチンはcall文を使って呼び出します。しかもそのとき引数を使って
変数の受け渡しをします。
たとえば
CALL RATE(B,TBASE,TAX,A,T)
で
SUBROUTINE RATE (E, TBASE, ANS, B, R)
に飛びます。ここでCALL RATE(B,TBASE,TAX,A,T)のBは,所得税の税率区分の数つまり５段階、TBASEには課税所得が、TAXには所得税額が、Aには課税所得の区分を示す配列変数,0,330,900,1800,3000,1000000000が、Tには税率0.1,0.2,0.3,0.4,0.5が入ります。これらの変数はメインプログラムからCALLで呼び出されたサブルーチン（副プログラム）へ引き渡されます。ここで注意して欲しいのはCALL文の中の引数とSUBROUTINE文の中の引数の名前が必ずしも一致していないことです。たとえばCALL文ではBですがSUBROUTINE文ではEですね。もちろん引数の名前は両者とも一致させても動きます。でも必ず一致させないといけないのであれば、このサブルーチンだけを部品として使った場合、メインプログラムとサブルーチンの間で変数名を統一させねばなりません。これは非常に厄介な作業です。そこでFOTRANでは、引数は変数名で対応するのではなく、単に順番のみで対応させています。つまりCALL文の１番目の引数がSUBROUTINE文の１番目の引数に対応し、CALL文の２番目の引数がSUBROUTINE文の２番目の引数に対応するわけです。
専門用語でいえば、グローバル変数とローカル変数の使い分けが可能になっているのです。かりにサブルーチンのなかで、メインプルグラムの中で使用している変数名と同じものを使用しても、メインプログラムの計算には何ら影響はありません。たとえば
real A,B,C
A=10
B=20
CALL TASIZAN(A,B,C)
WRITE(*,*) A,B,C
END
SUBROUTINE TASIZAN(E,F,G)
REAL E,F,G,C
G=E+F
C=100
RETURN
END
のプログラムだとサブルーチンの中でC=100としてますがこれはメインプログラムでの計算結果の表示に何ら影響を与えません。もし引数の変数名を同じにしてしまうと
real A,B,C
A=10
B=20
CALL TASIZAN(A,B,C)
WRITE(*,*) A,B,C
END
SUBROUTINE TASIZAN(A,B,C)
REAL E,F,G,C
C=A+B
C=100
RETURN
END
だと、サブルーチンでのc=100が意味をもってきます。各自確かめてね。この引数の使い方によってサブルーチンは汎用性の高い部品になります。所得税の計算プログラムの中でRATEというサブルーチンが３カ所から呼び出されているの気がつきましたか？まず
CALL RATE(E,WAGE,WAGEDE,WB,WR)
で呼び出されます。これは給与所得控除を計算するためです。給与所得控除は年収1000万円なら
180×0.4+(360-180)×0.3+(660-360)×0.2+(1000-660）＝220万円ですね。この計算方法は実はIF文を覚えようの例題で使った課税所得から税額を求めるときと似ていることに気づきませんか？次に
CALL RATE(S,WAGE,PREMI,PB,PR)
で再びRATEというサブルーチンが呼び出されます。今度は社会保険料控除の計算に使います。社会保険料は収入によって異なるのですが、大蔵省は課税最低限の国際比較などに使用するため便利な簡易計算方式を開発しています。これはあくまでも近似値として使われるもので、実際の納税の際には使われませんので念のため。大蔵省によると社会保険料は500万円以下の年収の人は年収×７％、５００万円を超える年収の人は500万円以下の部分は７％、500万円を超えた部分に２％をかけるとほぼあてはまりがいいそうです。私はシミュレーション分析にはいつもこれを使ってます。最後に
CALL RATE(B,TBASE,TAX,A,T)
で３度呼び出されます。今度は課税所得に対して所得税額を計算するのに使われます。
つまり、このRATEというサブルーチンは、所得税の計算のように段階的に計算していくときに便利な部品になっていることがわかりますか？
では最後に問題です。この所得税額を計算するプログラム自体をサブルーチン化して、年収５００万円から１０００万円まで１００万円刻みで増加させて、それぞれの所得税額を計算するプログラムを作成してみてください。
コンピュータによるシミュレーションで最も便利なものが収束計算です。収束計算とは、簡単にいうと解を見つけるのに、数式を解くのではなくて、片っ端から計算していけばそのうち正解にたどり着くでしょうという極めて粗雑なやり方です。よくミクロ経済学に出てくる制約条件付き最大化問題は、通常ラングランジュ乗数法を用いてかっこよく解きますが、実は、制約条件をみたす数値の組み合わせを全部試せば効用最大化の点は簡単にみつかります。ただし、人力でおこなうとやたらと時間がかかります。でもコンピュータは、無数に近い組み合わせを一瞬のうちにいやがりもせず計算してくれます。そこで例をあげてプログラムを組んでみましょうか・・・
まず、次のような問題を考えます。
効用関数 U=ＸＹ
消費者の予算制約 1000=100Ｘ＋100Ｙ
要するに1000円もってみかん（X財）とりんご（Ｙ財）を買いにいきました。効用を最大にするにはみかんとりんごをいくつ買ってかえればいいでしょう・・・というミクロ経済学の初歩でならう問題です。ラグランジュ乗数法を使えばみかん５個、りんご５個ということがわかります（この程度ならラグランジュ乗数法を使うまでもないですが・・・）。
この問題をフォートランプログラムで解いてみましょう。
REAL X,Y,U,P1,P2,M,ST
INTEGER I
P1=100
P2=100
M=1000
U=0
ST=1
X=1
10 Y=M/P2-(P1/P2)*X ! 下のGOTO文との間で繰り返し計算がおこなわれます。
U=X*Y !
WRITE(*,*) X,Y,U !
X=X+ST ! 数学的には変な式ですが、代入文です。Xの値が増えていくのがわかりますか？
IF (X.GT.9) STOP ! この式でｘ＞９を超えたときに繰り返しを止めます。
GOTO 10 ！10行目に戻るという命令です
STOP
END
このプロガラムを走らせると
1.000000 9.000000 9.000000
2.000000 8.000000 16.000000
3.000000 7.000000 21.000000
4.000000 6.000000 24.000000
5.000000 5.000000 25.000000
6.000000 4.000000 24.000000
7.000000 3.000000 21.000000
8.000000 2.000000 16.000000
9.000000 1.000000 9.000000
となります。３列目の数字が効用の値です。りんごとみかんをそれぞれ５個づつ買えば、効用の値は２５となり最大値になってますよね。ここが答えです。でもこれだと計算結果をすべて表示しないと答えがわかりません。そこで次のプログラムをみてください。
DOUBLE PRECISION X,Y,U,P1,P2,M,ST,UMAX
INTEGER I
5 FORMAT(F8.4,F8.4,F8.4,F8.4,F8.4)
P1=100
P2=100
M=1000
UMAX=-1E+10
U=0
ST=1
X=1
10 Y=M/P2-(P1/P2)*X
U=X*Y
IF (U.GT.UMAX) UMAX=U ！ここがポイント
WRITE(*,5) X,Y,U,UMAX,ST
IF (U.LT.UMAX) THEN ！ ここもポイントです。
X=X-ST ！ 何をしているかわかりますか？
ST=ST/10 ！
IF (ST.LT.0.00001) STOP ！
X=X+ST ！
GOTO 10 ！
END IF
X=X+ST
IF (X.GT.10) STOP
GOTO 10
STOP
END
このプログラムを走らせると
1.0000 9.0000 9.0000 9.0000 1.0000
2.0000 8.0000 16.0000 16.0000 1.0000
3.0000 7.0000 21.0000 21.0000 1.0000
4.0000 6.0000 24.0000 24.0000 1.0000
5.0000 5.0000 25.0000 25.0000 1.0000
6.0000 4.0000 24.0000 25.0000 1.0000
5.1000 4.9000 24.9900 25.0000 0.1000
5.0100 4.9900 24.9999 25.0000 0.0100
5.0010 4.9990 25.0000 25.0000 0.0010
5.0001 4.9999 25.0000 25.0000 0.0001
5.0000 5.0000 25.0000 25.0000 0.0000
となります。ここでわかるようにフォートランで求めた解は近似値です。たとえば小数点以下３けたまでもとめるなら、
最後から２行目の5.0001 4.9999 の組み合わせで止まります。それでは問題です。フォートランプログラムで住宅ローンの返済の計画を立ててください。あなたは1000万円を銀行から借りたとしましょう。金利が５％、返済期間が３０年だとすると、毎年の返済額をいくらにすれば借金が返せますか？毎年の返済額は同一になるようにしてください。
実践編１で連載を終了したつもりだったのですが、肝心なことを忘れてました。実際にシミュレーションプログラムを作成するには欠かせないデータの入出力です。合計を出してみようで配列変数にデータをセットするときにdata文を使用しましたね。データ数が少ない時はこれで十分なのですが、大量のデータをセットするにはあまり向いていません。そんなときはあらかじめ外部ファイルとしてテキスト形式のデータを用意しておいて、読み込ませる方が便利です。ほとんどの人がデータ処理をするにはエクセルのような表計算ソフトを使用していると思います。そこで、フォートランで読み込み可能な形式のテキストファイルをにエクセルで作成する方法を伝授しましょう。まずエクセルであらかじめデータを用意しておきましょう。たとえばこんな感じのファイル（エクセルで見てね）をね。全部で３８個のデータです。これは平成１０年の家計調査の年齢階級別の世帯主勤め先収入のデータを１歳刻みに加工したものです。このエクセル形式のファイルのままでは、フォートランのプログラムでは読み込めません。そこでエクセルでこのファイルをテキスト形式で保存する必要があります。ただし、ここで注意して欲しいのが、通常はデータ系列以外にタイトルとかデータの年数を示す文字列が入力されているはずです。でもこれらの文字列はフォートランに引き渡すとあとで面倒なので必ずこのファイルのように数字だけのファイルにしておいてください。
数字だけのファイルは、エクセルで、ファイル→名前を付けて保存を選択し、ファイルの種類の入力欄右の下向きの矢印部分をマウスでクリックして、テキスト（タブ区切り）.txtを選択してください。csv（カンマ区切り）.csvでもかまいません。今回はテキスト形式を使いましょう。ここではwage.txtというファイルを作成してみました。これでデータの準備はできました。あとはフォートランで読み込むだけです。ただし、このデータはフォートランのプログラムと同じディレクトリに保存しておいてくださいね。別の場所においておくことも可能ですが、ややこしいのでお勧め出来ません。
REAL WAGE(1:38) !
wageという配列変数に３８個のデータをセットすることを宣言
INTEGER I,OS !
整数型の変数としてIとOSを使用することを宣言
OPEN(UNIT=1,IOSTAT=OS,FILE='wage.txt',STATUS='OLD')
!wage.txtというファイルを装置番号１番として使用することを宣言
READ(1,*) (wage(I),I=1,38)
!装置番号１番のデータを読み込み
close(unit=1)
DO I=1,38
WRITE(*,*) I,WAGE(I)
END DO
STOP
END
このプログラムにおいて気をつけることは使用する変数の宣言を忘れないことです。まず、読み込んだデータをフォートランで使用するために配列変数を宣言します。ここではREAL
WAGE(1:38)においてWAGEという名前の配列変数を添え字１から３８までを実数型の変数として使用することを宣言しています。ここまでは常識ですが、忘れやすいの次の宣言文
INTEGER I,OS
の部分です。これは整数型の変数としてIとOSを使用するという宣言文です。Iという変数は、配列変数の添え字を意味する変数として使用します。配列変数の添え字はかならず整数型で宣言してください。実数型でも問題なさそうですが、コンピュータの数値計算においては、整数（つまり小数の出てこない数字ね・・・）の方が計算スピードが速く、メモリ（記憶容量）も消費しないからです。メモリの消費は最近のパソコンではほとんど無視できますが昔のパソコンでは重要だったのですよ。なにしろわたしの最初に購入したパソコンは６４キロバイトのメモリ空間しか使用出来なかったぐらいです。次のOSという変数は、外部ファイルを読み込むときに使用する命令文であるOPEN文のなかで使用します。IOSTAT=OSの部分です。これの意味は知りません。決まりごとですから・・・。
OPEN(UNIT=1,IOSTAT=OS,FILE='wage.txt',STATUS='OLD')
ですが、もし読み込みたいファイルが複数あった場合には
OPEN(UNIT=1,IOSTAT=OS,FILE='wage.txt',STATUS='OLD')
OPEN(UNIT=2,IOSTAT=OS,FILE='wage2.txt',STATUS='OLD')
というように、UNITの数をかえて複数記述してください。ただし、UNIT番号は５番までにしてね。なぜか私の使用しているフォートラン９０では６番目からはエラーがでることがあります。フォートラン７７のときは１０個ぐらい指定しても大丈夫だったのですが。これがフォートラン９０固有の制限なのかそれともdec版の固有の制限なのかは知りません。と書いていましたがやっと６個以上のファイルを一度に指定してもエラーがでない方法がわかりました。わたしがこれまでに作成していたプログラムでは最大６０個のファイルを一度に指定しているのですが、そのときにはエラーが出ませんでした。実はclose文の場所が大事だということがわかりました。read文とclose文の間にwrite文が入っているとopen文で指定しているファイルを破壊してしまいます。だから、close文はread文の直後に記述する必要があったのです・・・
もし同時に６個以上のデータ系列を読み込みたいときは、一旦配列変数への読み込みが終わった段階で、UNITをクローズします。外部ファイル読み込みのために確保した領域を白紙に戻してやるわけです。たとえば
DO I=1,5
CLOSE (UNIT=I)
END DO
を記述すれば５個のUNITをまとめて白紙に戻せます。STATUS='OLD'の部分はすでに作成済みのファイルの読み込みをおこなうことを意味しています。古いファイルですね。データの読み込みは
READ(1,*) (wage(I),I=1,38)
の部分です。もしUNIT=2のファイルを読み込むときは、READ(2,*)にします。この
(wage(I),I=1,38)
の部分は（）でくくるだけでDO文を使わずにループ（繰り返し）を実現しています。通常のループ文を使用して
DO I=1,38
READ(1,*) wage(I)
END DO
と記述してもかまわないのですが、この方が１行で記述できるので見やすくなってます。
DO I=1,38
WRITE(*,*) I,WAGE(I)
END DO
の部分は外部ファイルとして読み込んだファイルを画面に表示させているだけです。
今回は、フォートランプログラムで計算した計算結果を活用する方法を考えてみましょう。これまでの連載では計算結果は、dos窓で表示させていました。これでもマウスでカットアンドペーストすればエクセルなど表計算ソフトに取り込むことができます。でももう少しスマートなやり方があります。それは、テキスト形式のファイルに出力することです。前回、テキスト形式のファイルの読み込みをしましたね。今度はその逆をおこないます。
REAL WAGE(23:60),tax(23:60)
INTEGER I,OS! 整数型の変数としてIとOSを使用することを宣言
OPEN(UNIT=1,IOSTAT=OS,FILE='wage.txt',STATUS='OLD')
!wage.txtというファイルを装置番号１番として使用することを宣言
READ(1,*) (wage(I),I=23,60) !装置番号１番のデータを読み込み
CLOSE (UNIT=1)
DO I=23,60
tax(I)=0.1*WAGE(I)
END DO
OPEN(UNIT=2,IOSTAT=OS,FILE='tax.csv',STATUS='new')
do i=23,60
write(2,*) I,',',tax(i)
end do
CLOSE (UNIT=2)
STOP
END
前回のプログラムに少し手を加えたものですが、配列変数の添え字の下限と上限もついでに変更しておきました。
WAGE(23:60)
としていますね。これは、wage.txtのデータは、２３歳から６０歳までの給与収入のデータなので、下限を２３、上限を６０にしてみました。この方が直感的にわかりやすいと思います。つまり配列変数の添え字は別に途中からつけてもかまいません。今回はこの２３歳から６０歳の給与収入に税率１０％の比例税を課したときの税額を計算して、計算結果をファイルに出力するプログラムを作成しています。
DO I=23,60
tax(I)=0.1*WAGE(I)
END DO
の部分が給与収入に１０％の比例税をかけて、税額をtaxという変数に代入しているところです。このtaxというデータをファイルに出力するには、あらかじめファイル出力の宣言が必要になります。
OPEN(UNIT=2,IOSTAT=OS,FILE='tax.csv',STATUS='new')
では、tax.csvという名前のテキスト形式のファイルを出力することを宣言しています。これは、テキスト形式のファイルを読み込む場合の宣言と似てますね。でもよくみるとSTATUS='new'のところが違います。読み込みのときは、STATUS='old'でした。つまり古い既存のファイルでなく、新しい新規のファイルを宣言しているわけです。それともうひとつ注目して欲しいのファイル名の拡張子を
.csvにしています。通常テキストファイルを意味するファイルの拡張子には.txtを使用します。.csvという拡張子は実は、カンマで区切られたテキスト形式のデータを意味しています。カンマで区切ったテキストデータは、エクセルなどの表計算ソフトで簡単に読み込むことができる便利な形式です。
do i=23,60
write(2,*) I,',',tax(i)
end do
が実際に、カンマで区切られたデータを書き出しているところです。','を使ってカンマをデータIとtax(i)の間に挿入しています。
Copyright(c) 1999-2001 by Kyoji Hashimoto
2001/04/17 19:46:23