Pythonを使った経済分析について

1
3
P
y
t
h
o
nを使った経済分析について
菅
準一
MATLABや P
y
t
h
o
nを使った経済学動学のテキストを手にしたとき、経済学教育
において数値計算やプログラミングを教える際の障害となってきたもののいくつか
巴m
a
t
i
c
aや MATLAB を使った演
を取り除けるのではと感じた。障害の一つは Math
習だと、卒業後は自分でこれらのソフトを購入しないと使い続けることができない。
個人で高額なソフトを購入することまで考えると、 M
a
t
h
e
m
a
t
i
c
aや MATLAB に依
存したスキルを積み上げていくことに二の足を踏む。二つ目の障害は使い勝手は我
慢することにして、クローンの Maximaや O
c
t
a
v
eを使えば良いが、経済学におい
て数式処理や数値計算をすること、つまりコンピュータを使って数学すること
(
d
o
i
n
gmath
巴m
a
t
i
c
sb
yc
o
m
p
u
t
e
r
)に重点をおいてスキルを磨いていったとき、興味が
他の分野に移ったときの対応を考える必要がある
O
同じ無料のソフト(こういう言
Lミ方は適切でないことは承知しているが、問題を際立たせるレトリックだと考えて
y
t
h
o
nは、様々なモジュールが用意されており、コンビュータを使って
欲しい。) P
数学すること (
d
o
i
n
gm
a
t
h
e
m
a
t
i
c
sbyc
o
m
p
u
t
e
r
)にとどまらない広い応用が考えられ
る
。
キーワード:Python，EconomicDynamics，Linux，Sag巴:OpenS
o
u
r
c
eMathematics
S
o
f
t
w
a
r
巴
日次
1 はじめに
2 町
北honを対話モード (
i
n
t
e
r
a
c
t
i
v
emode)で使う
3 同社10nの NameE汀 orについて
4 町
北honでの繰り返しブロック:インデント
5 p)ぇhon のスクリプトファイ)~ s
i
n
x
_
u
t
L
8
.
p
yの作成と実行
6 乃
花honを使った経済分析
7 Sage:OpenSourceMathematicsS
o
f
t
w
a
r
eを使う
8 おわりに
AppendixL
inuxT
i
p
s
Vo.
l10
1
4
NO.2
1 はじめに
殴米では認知度の高い言語である MATLABや Pythonを使った経済学動学のテキスト 1を
子にしたとき、経済学教育において数値計算やプログラミングを教える際の障害となってき
たもののいくつかを取り除けるのではと感じた。障害の」つは Math巴maticaや MATLABを
使った演習だと、卒業後は自分でこれらのソフトを購入しないと使い続けることができな Lし
個人で高額なソフトを購入することまで考えると、 Mathematicaや MATLAB に依存したス
キルを積み上げていくことに二の足を踏む。それではクローンである Maximaや Octaveを
使えば良いが、ソフトをダウンロードしてインストールした経験がないものにとっては抵抗
があるようだ。苦労してインストールして使ってみたが、本家本元のほうが何かと初学者に
は使いやすい。二つ日の障害は使い勝手は我慢することにして、クローンの Maximaや
Octave を使えば良いが、経済学において数式処理や数値計算をすること、つまりコンビュ
ータを使って数学すること (doingmathematicsbycomputer)に重点をおいてスキルを磨いてい
ったとき、興味が他の分野に移ったときの対応を考える必要がある。同じ無料のソフト(こ
ういう言い方は適切でないことは承知しているが、問題を際立たせるレトリックだと考えて
欲しい。) Pythonは、様々なモジュールが用意されており、コンピュータを使って数学する
こと (doingmathematicsbycomputer)にとどまらない広 LリJ吾用が考えられる。 2 ここでは第 1
に Python という言語の使いやすさ、習得の容易さを簡単な例を用いて検討した
O
第 2に
EconomicDynamicsのホームページにある Pythonのコードと MATLABのコードとを比較し
ながら、 Pythonを経済学教育に使う可能性についても検討した。
2
Pythonを対話モード(in
t
e
r
a
c
t
i
v
emode)で使う
新しい電気製品を買ったとき、時には類似の製品から類推して使えることもあるが、取
り扱い説明書を片子に子順を踏んで使うのが最もトラブルが少ない。しかし何となく使えそ
うな気がして、ボタンを押してしまうのも人の性かもしれない。ターミナルから python と
入力して起動すると、 pythonの入力プロンプトである
>>>(chevron というようだ) (こ変わ
る。通常は四則演算を試したりするのだが、データを可視化することが理解を深めることが
多いので、 Pythonのグラフィック機能がどのようなものか試したい。まずサインカーブを
-2π から 2π の範囲で描かせようと思う
o
Mathematica流に Plot[Sin[
x
l，{
x，2句 i，
2*Pi}l
なのか、それとも MATLAB 流に x=-2*pi:pi/20:2*pi; plot(x，sin(x)) なのか。
Mathematicaの関数の書き方はかなり特殊なので、 MATLAB流を試す。ところが:で
J
o
h
nS
t
a
c
h
u
r
s
k
i，
2
0
0
9，
E
c
o
n
o
m
i
cD
y
n
a
m
i
c
s，
MITP
r
e
s
s
例えば ZOPE
，I
n
k
s
c
a
p
eなど
P
y
t
h
o
nを使った経済分析について
S
y
n
t
a
x
E
r
r
o
rが出てしまう
O
1
5
試しに p
i と打つと NameError:name'
p
i
'i
sn
o
td
e
f
i
n
e
d となり、
s
i
n
(
x
) と打つと NameError:nam
巴'
s
i
n
'i
sn
o
td
e
f
i
n
e
d となる O 特別なことをしなくても、 p
iや
smが使える Mathematicaや MATLABのようなソフトとは違うようだ。まずはマニュアルを
読んだ方が無難なようだ。数学でよく使う定数 π や関数 s
i
nx を使うためには、 math とか
numpyというモジュールが必要らし t
¥
' グラフを描くには、 y=s
i
nx という関係を満たす点
0
(
x，
y
) の集まりを用意しないといけないが、まず司2
π から 2π の範囲で x座標の点のリスト
を作るには arangeも必要になるようだ。さらに p
l
o
tさせるには、 p
l
o
t用のモジュールが必
要であることが分かる。 m
a
t
p
l
o
t
l
i
bの p
y
l
a
b というモジュールが高機能であることが分かっ
た。まず必要なモジュールを i
m
p
o
r
tコマンドを使ってインポートする。
>>>from numpy import pi，sin，
arange
show
>>>from py1ab import p1ot，
もし次のようなエラーメッセージが出るようだと、適当なサイトからファイルをダウンロ
ードしてインストールする作業が必要となる O
Traceback (most recent ca11 1ast):
Fi1e "cstdin>"， 1ine 1， i
n cmodu1e>
ImportError: No modu1e named py1ab
a
t
p
l
o
t
l
i
bのサ
p
y
l
a
bというモジュールがコンビュータにインストールされていないので、 m
イトで、自分が使っている
o
s用のバイナリーファイルをさがしてモジュールをインストー jし
する。その際使っている pythonのパージョンに対応したものであることを確認する必要が
ある。
必要なモジュールの importがうまくいったなら、まず、 a
r
r
a
y
r
a
n
g
e
(
s
t
a
r
t
=
O，
s
t
o
p，
step=l，
t
y
p巴
cod
巴
=Non
りあるいは短縮形の a
r
a
n
g
巴を使ってグラフを描かせるための x軸上の点をつ
2
π から 2π の区間を L 刻みに細分する。
くる。 20
>>>x=arange( 2*pi，2*pi，
pi/20)
戸
とすると、 2π から 2π をよL刻みで等分した 8
0個の点が得られる。 3
20
3 似た関数に、
r
a
n
g
eがあるが、整数しかとれないので注意
Vol
.1
0
1
6
NO.2
>>>plot(x，
sin(x))
>>>show()
とすると、先の 8
0個の点 x について sin(x) が計賀されて、別ウインドウが立ち上がり、サ
インカーブが描かれる
O
しかし、 2
π までちゃんとグラブが描かれていないことに気付く。
>>>x
array([-6.28318531，
6.12610567])
仕方がないので
>>>x=arange(-2*pi，
2*pi+pi/20，
pi/20)
array([-6.283l8531， ••. ，
6
.2831853l])
となり所望の結果が得られる。
まとめ
Pythonをターミナルから起動する (pythonでも i
d
l
巴
， i
pythonなどと打って好みの環境を起動
するとよ P)、下のリストを試したいとき、入力プロンプト>>>より後の部分が入力する部
分である。>>>を入力する必要はない。それじゃ一入力するぞと思って、気持ちの区切りと
ともに、うっかりスペースを打ったりするとエラーになるので注意が必要だ。 Pythonにと
ってスペースを打つことつまり ind叩 tすることは、大きな意味を持っているからだ。
>>>from numpy import sin，
pi，
arange
>>>from pylab import plot，
show
>>>x=arange(ー2*pi，
2*pi+pi/20，
pi/20)
sin(x))
>>>plot(x，
>>>show(
)
実質 2行ですみます。何といっても、グラフを描くために x軸上に取った点のリスト(正確
に言うと配列 (array)) に sinOが適用できることが強みです。 4
高価な Mathematicaで Plot[Sin[x]ム -2*Pi，
2叩 i
] でサインカーブを描かせたものと、 Python
P
y
t
h
o
nを使った経済分析について
1
7
のものとを比較すると以下のようになる。左側が Mathematicaのものである
O
ソフトウェア
d
o
i
n
g
任せの安直さという点では Mathematicaの方が、コンビュータを使って数学をする (
m
a
t
h
e
m
a
t
i
c
s
)ことを始めるには良いかも知れない。
LO
05
1
.
0
3 Pythonの NameErrorについて
i と入力すると
すこし寄り道になるが、 NameErrorが気になる。 python を起動直後に p
NameError:name'
p
i
'i
sn
o
td
e
f
i
n
e
dとなるの fromnumpyi
m
p
o
r
ts
i
n，
p
iとした後ではがを使って
もエラーが出ない。少し調べてみると、現在使える名前 (Name)のリストを列挙するには
d
i
r
Oでよいらし Lミ。そこで py出 on起動直後に d
i
r
Oと打ってみると
[
' builtins
doc
一
一name
package '
]
確かに p
iや smはない。そこで fromnumpyi
m
p
o
r
ts
i
n，
p
i とした後で d
i
r
Oで確認すると、
，
[
builtins
doc
name
package
'
p
i
'，，
s
i
n
'
] となってお
iや smという名前が使える。
り
、 p
ネヅトで調べていて l
m
p
r
o
tnumpyという方法があってそのときは n
u
m
p
y
.
p
iとか numpy.sm
としないと使えないという記述を見たことがある。 d
i
r
Oで、探ってみよう。
p
y
t
h
o
nを再起動して i
m
p
o
r
tnumpyの後 d
i
r
Oとすると、
，
[
builtins
doc
目
name
package
'
n
u
m
p
y
'
] となっており、 p
i
や smという名前は使えない。しかし以下で示すように n
u
m
p
y
.
p
iや n
u
m
p
y
.
s
i
nは使える。
>>> n
umpy.pl
.
4 リスト(l
i
s
t
)は違う
t
y
p
eのデータをまとめることができるけれど、自己列 (
a
π
a
y
)は同じ t
y
p
eのものしかまとめ
られないという制限がある c
Vo
l
.1
0
1
8
No.2
3.141592653589793
>>> numpy.sin{numpy.pi/2J
1
.0
numpy.plや numpy.sm はめんどくさいので、かわりに p
i
.
s
i
n を使いたければ、以下で示すよ
うにそのように定義すればよい。
>>> pi=n吐血py.pi
>>> pl
.
3.141592653589793
>>> sl
.n
=numpy.s工n
>>> sin{pi/2)
1
.0
>>> dir(
)
ー
[
builtins
doc
name
'__package__'， 'numpy'， '
p
i
'， 'sin']
また、 fromnumpyimport*どしたとき何がおきるか dirOで調べてみるどいい。
少し話を戻しますと、 fromnumpyimportや frompylabimportでエラーがでなかったのは、
あなたかコンビュータの管理者が numpyと matplotlibをダウンロードしてインストールして
いたからです。 Pythonの標準インストールではこの 2つのモジュール(ノ tッケージ)は組
み込まれていません。いま起動している Pythonの標準組み込みモジュールはどんなものか
知りたいときには、以下のようにすればよ p 。
>>>import sys
>>>sys.builtin module names
{
' builtin
main
' ast'， 'bisect'， 'codecs'，
codecs cn'，
，codecs hk'， 'codecs iso2022'， 'codecs jp'， 'codecs kr'， 'codecs tw'，
'collections'， 'csv'， 'functools'， 'heapq'， 'hotshot'， 'io'， 'json'，
，locale'， '_lsprof'， '
_
md5" '
_
multibytecodec'， 'random'， 'sha'，
，sha256'， 'sha512'， 'sre'，
I
warningsI I
I
struct'，' subprocess'， 'symtable'，
weakrefI ， I winregI ， Iarray1 I
audioopI ， Ibinascii1 ，
I
'cPickle'， 'cStringIO'， 'cmath'， 'datetime'， 'errno'， 'exceptions'，
'future builtins'， '
g
c
'， 'imageop'， 'imp'， 'itertools'， 'marshal'，
P
y
t
h
o
nを使った経済分析について
嘗
1
9
math'， '
m
m
a
p
'， 'msvcrt'， '
n
t
'， 'operator'， 'parser'， 'signal'， '
s
t
r
o
p
'，
'
s
y
s
'， 'thread'， '
t
i
m
e
'， 'xxsubtype'， 'zipimport'， 'zlib')
a
t
p
l
o
t
l
i
bは標準組み込みモジュールとして名前があがっていないので自分でイン
numpyや m
ストールする必要がある。それに対して s
y
sはいま起動している Pythonの標準組み込みモ
ジュールとして、名前があがっている。したがって外部からインストールしなくても import
s
y
sですぐに使い始めることができる
O
標準組み込みモジュール s
y
s とそれ以外のモジュール numpy，
p
y
l
a
b との遣いは理解できた
と思うが、 d
i
r
Oて何なのか気になる。 i
m
p
o
r
t しないで直ぐに使い始めた。実は d
i
rは標準組
み込み関数である。標準の組み込み関数や変数は、 dirO としただけでは見られないように
なっている
O
>>>
dir(
)
，
[
builtins
doc
name
実は d
i
r
(
)は上の出力の先頭にある
>>>
builtins
package
'numpy'， '
p
i
'， '
s
i
n
'
l
の中にある。実際
dir( builtins
['
A
r
ithmeticError'， 'AssertionError'， '
A
ttributeError'， 'BaseException'，
'BufferError‘ 'BytesWarning'， 'DeprecationWarning'， 'EOFError'， 'Ellipsis'，
'EnvironmentError'， 'Exception'， 'False'， 'FloatingPointError'，
'
F
utureWarning'， 'GeneratorExit
'， 'IOError'， 'ImportError'， 'ImportWarning'，
'IndentationError'， 'IndexError'， 'KeyError'， 'Keyboardlnterrupt'，
'LookupError'， 'MemoryError'， 'NameError'， '
N
o
n
e
'， 'Not工mplemented'，
'NotlmplementedError'， 'OSError'， 'OverflowError'，
'PendingDeprecationWarning'， 'ReferenceError'， 'RuntimeError'，
'RuntimeWarning'， 'StandardError'， 'Stoplteration'， 'SyntaxError'，
'SyntaxWarning'， 'SystemError'， 'SystemExit'， 'TabError'， 'True'， 'TypeError'，
'UnboundLocalError'， 'UnicodeDecodeError'， 'UnicodeEncodeError'，
'UnicodeError'， 'UnicodeTranslateError'， 'UnicodeWarning'， 'UserWarning'，
'ValueError'， 'Warning'， 'WindowsError'，
'ZeroDivisionError'，
debug
doc
import
name
a
b
s
'， '
a
l
l
'， '
a
n
y
'， 'apply'， 'basestring'， '
b
i
n
'， '
b
o
o
l
'，
'__package__'， '
Vo.
l1
0
2
0
NO.2
'buffer'， 'bytearray'， 'bytes'， 'caユlable'， 'chr'， 'classmethod'， 'cmp'，
'coerce'， 'compile'， 'complex'， 'copyright'， 'credits'， 'delattr'， 'dict'，
'dir'， 'divmod'， 'enumerate'， 'eval'， 'execfile'， 'exit'， 'file'， 'filter'，
'float'， 'format'， 'frozenset'， 'getattr'， 'globals'， 'hasattr'， 'hash'，
'help'， 'hex'， 'id1 ， 'input'， lint'， 'intern'， 'isinstance'， 'issubclass'，
'iter'， '1en'， 'license'l 'list'， 'locals'， 'long'， 'map'， 'max'， 'memoryview'，
'min'， 'next'， 'object'， 'oct'， 'open'，目。 rd'， 'pow'， 'print'， 'property'，
'quit'， 'range'， 'raw_input'， 'reduce'， 'reload'， 'repr'， 'reversed'I 'round'I
'set'， 'setattr'， 'slice'， 'sorted'， 'staticmethod'， 'str'， 'sum'， 'super'，
tuple'， 'type'，
・ unichrI ， unicode ， varsI ， xrangeI ， zip
，
I
d
i
rが
>>>
builtins
1
I
I
I
1]
中に含まれていることが分かる。
builtins .dir()
[
' builtins
doc
name
package
p
i
'， 'sin']
'numpy'， '
色々ややこしい話をしたが、ここでいう名前とは~、下に示すように varsO で表示される辞書
(
d
i
c
t
i
o
n
a
r
y
)のキー (
k
e
y
)の部分であり、その実態は:の次にある値 (
v
a
l
u巴)なのだ。
>>> vars()
{
' builtins
目
一
一
一package
・ <module'
builtin ' (built-in)>，
'
: None，
name '
:
main
'pi': 3.14ユ592653589793，
'numpy': <module 'numpy' from '/Library/Frameworks/.../ init .pyc'>，
doc '
: None，
sin1 :
I
< 吋u
nc
builtins
-}
l
の
'sin'>}
の値は <module
builtin
builtin ーであり、 p
lはすでに
かる。さらには、
(built-in)>つまり標準組み込みモジュ
3
.
1
4
1
5
9
2
6
5
3
5
8
9
7
9
3という{直を持っていることが分
S
l
l
l は数値だけだなくリストにも適用でき、
u
f
u
n
c
(
u
n
i
v
e
r
s
a
lf
u
n
c
t
i
o
n
)である
ことも確認できる c スクリフトを書いて実行させているとき、時々 dirOでどのような名前
が使えるのか確認するといい。
Pythonを使った経済分析について
4
21
Pythonでの繰り返しブロック:インデン卜
繰り返しというと、 f
o
rを思いつく人も多いと思う。繰り返しの回数を計算式やその他の
形で指定できる場合は使い勝手がよいが、どれくらいの繰り返し回数になるか計算が可能だ
けれど厄介なときとか、そもそも何回繰り返せばよいか前もって分からない場合などには、
whileを使う子がある。 5
横軸の点と縦軸の点をおさめる x
，
yという空のリスト(オブジェクト)を生成する。
>>>x=[1
>>>y=[ 1
たとえば xのタイプを typ巴コマンドで確認するには
>>> type{x)
<type 'list'>
さらに d
i
r
(
x
)とすると
>>> dir{x)
，
[
add
c1ass
de1s1ice
doc
getattribute
iadd
mu1
repr
sets1ice
contains
de1attr
de1item
eq
format
ge
一
一 getitem__
gets1ice
gt
imu1
ne
init
new
reversed
sizeof
iter
reduce
rmu1
str
hash
1e
1en
1t
reduce ex
setattr
setitem
subc1asshook
'append'，
'count'， 'extend'， 'index'， 'insert'， 'pop'， 'remove'， 'reverse'， 'sort']
u
i
l
t
i
nメソッドが表示される O
append以下にオブジェクト xが使える b
5f
o
rを使っても w
h
i
eと同じことができる。十分大きいと思える繰り返しの回数を指定しておいて i
fで中断
する方法である c 開始点を s=-2*pi刻みを d=pi/20 と指定して、 x=[s+iサ f
o
rii
nxrang巴(
0，
10000)i
f
s
+
i
*
d
<
=
2
*
p
i
]というリストの内包表現を使う。 x
r
a
n
g
eは g
e
n
e
r
a
t
o
rで数 1 を一つ作っては f
o
rブロックの処理
をして、 i
+
1に置き換え処理を繰り返す c r
a
n
g
eだと l万個の整数からなるリストがつくられ、そのリスト
上を 1 が動いていくので、 x
r
a
n
g
eよりメモリを多く消費する
P
y
t
h
o
n
3では x
r
a
n
g
eは r
a
n
g
eに吸収されて
なくなる
l
.1
0
Vo
2
2
NO.2
以下のようにして点のリストを作る、オブジェクト x が使える出来合の append というメ
ソッドを使うと簡単である。このあたりは Pythonがオブジェクト指向の言語である強みで
あり、言語を習得することを容易にするかもしれない。神の後はコメントなので入力する必
要はない。ただ、これらをひとまとめにして pythonスクリプトとして別ファイルに保存す
るときスクリプトの可読性をあげるために必要なものもある。
>>>i=・2*pi
#区間の始まりを設定
>>>d=pi!20
非刻みを設定
>>>while i <=2*pi:
持2*pi を含むようにする。コロン:を忘れない
x.append(i)
軸プロンプトが...に変わる
y.append(sin(i))
while のブロックは半角 4文字分字下げを忘れない
非
i=i+d
時半角スペース 4つの後から入力
持プロンプトに対して何も入力しないで E
n
t
e
rするとブロックの終わり
プロンプトが>>>に戻ったらエラーがなく処理が終わったことになる。
そこで、
>>>plロt(x，
y
)
show(
)
5 Pythonのスクリプトファイル sinx_utC8.pyの作成と実行
対話モードで、使っているとき、その編集機能が使い慣れたエディタに及ばないことが多い。
特に日本語を使ったりするのは、使い慣れたエディタがいい。また対話モードで行ったこと
の全体を傭撤して改良を加えたりする場合には、これまで行ったコマンド入力のヒストリー
全体をエディタの画面で眺めてみるとよい o Python のデフォルトの IDLE で F
i
l
e
NewWindowで立ち上がるエディタでスクリプトを書くのが LミL功ミも知れない。これだと
Runメニューが表示され、修正しては再実行することができるので便利で、ある。最初はこれ
を使うといい。慣れてくると、対話的に使うための PythonS
h
e
l
lの winodw と使い慣れたエ
ディタ windowを並べておいて、最初の実行には PythonS
h
e
l
lから i
m
p
o
r
tを使って実行し、
h
e
l
lで r
e
l
o
a
dを使うとよい。スクリプト
エディタで修正し保存した後の再実行には PythonS
の可読性を高めるために、日本語のコメントを入れたファイル sinx_utL8.pyを utL8が読み
書き保存可能な設定を施したエディタで作成してみた。
Pythonを使った経済分析について
時 -2pi<=x<=2pi の範囲で sin x
2
3
のグラブを描く
from numpy import pi，
sin
#必要なものだけを読み込む
from pylab import plot，
show
持この後で dirOとするとよい
print(dir(
))
非import *との違いは
x=口
y=[
]
i=・2*pi
非区間の始まりを設定
d=pij20
持刻みを設定
while i<=2*pi:
持2*pi を含むようにする。コロン:を忘れない
hile のブロックは半角
4文字分字下げを忘れない
x.append(i)
村
y.append(sin(i))
持半角スペース 4つの後から入力
i=i+d
持ブロックの外だから字下げしない
print(x)
print(y)
plot(x，
y
)
非グラフを作図する
show()
韓これがないと画面表示されない
を作ってこれを importで読み込んで実行すると次のようなエラーメッセージがでた。
>>> import sinx uヒf 8
Traceback (most recent call last):
File "<stdin>"，line 1，in <module>
File sinx_utf_8.py"，line 1
四
SyntaxError: Non・且 SCII character I¥xe3・in file sinx_utf_8.py on line 1，
but no encoding declared; see http://www.python.org/peps/pep-0263.html
for details
非アスキー (Non-ASCII)文字が使われているが、エンコーディング方式が指定されていない
という指摘を受けた。コメントアウトしているのだから、何語で書こうと無視して欲しいと
思うがそうはいかないようだ。
そこで、文字コード方式を指定する行を第 1行に加えた。
非ー *
-coding:
utf 8 -*-
Vol
.1
0
24
NO.2
そして plot(x，
y
)の行の前に
print
1
1
日本語が表示できます"
を付け加えて、保存し実行する。
すると、エラーを出すことなく実行される。
以上まとめると、 Pythonのスクリプトを作成するためには、最初は IDLEのエディタを
使う、 Pythonの Keyword (予約語)を色付けして表示してくれたり、 forや whileなどのブ
ロックでは自動的に indentしてくれることや、エラーが出るとその箇所を表示してぐれるの
で便利である。日本語を使ったりして、足らない所が出てくると自分の使い慣れたエディタ
を使えるので、プログラムを開発する環境も整っている。
6 Pythonを使った経済分析
以下に EconomicDynamics からあまり難しくない問題を選んで Pythonのスクリプトを書
いてみた。
マルコフチェーンが t回推移した後に取る状態を確率変数 Xtとするとき、 Xt=yとなる
確率を求めるプログラムに、必要な数値を対話的に入力する部分を付け加えた。一応入力に
対して簡単なエラーチェックも付けてみた。 Python に慣れていない人が利用するために対
話的なプログラムを書いたものである。
持・*-
coding:utf 8 -*
ー
幽
神 Filename: testmarginal_input.py
持 Author: Junichi Suga
非 Date: August 2010
非 Corresponds t
o: Exercise 4.2.3
持frompylab import *
from mygenfinitemc import samp1e，MC 非 Import from listing 4.4
pH
((0.971， 0.029， 0.000)，
(0.145， 0.778， 0.077)，
(0.000，0.508， 0.492))
adic=，
[NG'，'MR'，'SR'
]
#辞書は勝手にソートされる
1，
2
1
bdic=[0，
非ソート前の状態を残す
Pythonを使った経済分析について
2
5
持辞書をつくる O 勝手にソートされる
state=dict(zip(adic，
bdic))
print 四マルコフチェーンで t回推移した後で状態 yが生じる確率を H
print "n同の実験で y が出た比率として計算するプログラムです。刊
print
1
1
1
1
print
print "初期状態の確率分布を入力します。"
c，
d
="11
，
0.0
while c!='y'or d!=l.O:
print("入力した
3つの数字の和が1.0 になるようにして下さ Lミ
C
11)
a=[
]
for i in range(3):
a.append(input("状態
"+adic[i]+" の確率はー>"))
d=sum(a)
print "psi="，
a
c=raw_input(祖入力した初期分布 psiは正しいですか (y/n)? "
)
持 Initial condition
psi=tuple(a)
t=O
while t<5 or 30<t:
print
t=i時叫 (
"
5から 30の範囲の整数を入力してくださ Lに推移の回数
t= ，
，
)
n=O
while n<10000:
n=input("サンプルの数
n>=10000，
n= "
)
"
y=闘
while y not in state:
神yが正しく入力されるかチェック
y=raw input(目確率を求める状態 y を半角英数大文字で入力して下さ Lに NG/MR/SR y= 目
)
h
MC(f=psi，
p=pH，X=sample(psi))
持
Create an instance of class MC
n，
state[y])
m=h.marginal(t，
" 回の推移の後で、状態が
print t，
になる確率は
，
I
l y，
，
" m，11
です。"
reload
以上を Pythonのスクリプトとして testmarginal input
.py に保存する。 testmarginal_input
.pyか
←
s として参照でき
ら必要なものを importする時に、名前が長くならないように、短縮形で t
2
6
Vol
.1
0
NO.2
るようにする。またエディタでファイルを手直しして、保存した後にも短縮した名前 t
sで
呼び出すことができる。 re!oad(ts) とすると修正が反映される。以下は実行した画面をコピ
ーしたものです。
>>>import testmarginal_input as ts
マルコフチェーンでヒ回推移した後で状態 yが生じる確率を
n 回の実験で y が出た比率として計算するプログラムです。
初期状態の確率分布を入力します。
入力した
3つの数字の和が1.0になるようにして下さい。
状態
NG
の確率は
ー
>0.3
状態
MR
の確率は
・
>0.4
状態
SR
の確率は
ー
>0.3
psi= [0.3，0.4，0.31
入力した初期分布 psiは正しいですか (y/n)? y
5から 30の範囲の整数を入力してくださ Lに推移の回数
サンプルの数
t= 24
n>=10000，
n= 10000
確率を求める状態 y を半角英数大文字で入力して下さ lに NG/MR/SR y= NG
24 回の推移の後で、状態が
NG
になる確率は
0.8041
です。
ここまでが、プログラムを対話的に実行した画面のコピーである
O
慣れてくると、この対話形式の入力が次第に煩わしくなる。ここまで対話的でなくても良
5回推移した後で状態が SRである確率を、標
いと思えるようになる。そうなってくると、 1
本数 5000で求めたくなったら、いちいち re!oad(同として対話的なプログラムを立ち上げる
よりも、オブジェクト指向プログラミング (OOP:ObjectOrientedProgramming) を活かして次
のようにするとよ t¥0
>>> ts.h.margina1(15，
5000，
ts.state['SR'1)
0.0336
まさに何をしているかは作った本人のみが知ることができる寡黙なプログラムである。乱数
Pythonを使った経済分析について
27
を発生させて計算しているから、結果は常に変動する O
つぎに単独で用いる場合にはその部分の処理が行われ、他のスクリプトから呼ばれる場合
は処理されないようにしたい場合は if
非
name ==' main ':ブロックを使う。
/usr/bin/env python
非・*・ coding:utf 8 ・
*
ー
幽
非 Joint Distributions
持 list on page 81
神 joint r utf
.py から呼び出して使うので、 import のときには
非 print させないように、 if 文を使っている。
持 from numpy import *
pH
((0.971，0.029， 0.OOO)，非デフォルトのカーネル
(0.145， 0.778， 0.077)，
(0.000， 0.508， 0.492)}
pSi
(0.2，0.2， 0.6)
非デフォルトの初期分布
state={'NG':
0，'MR':
1，'
SR':2} 非ステートの辞書 state を定義
X=('NG'，'MR'，'NG')
非デフォルトのパス
持確率を求めたいパスを指定するには
>>>path_prob(x=('NG'，'MR'，'MR'}}
psi=pSi，
x=X}:
def path prob(p=pH，
持キーで数値を呼び出す
prob=psi[state[x[Olll
for t in range(len(x}・1
}:
prob=prob*p[state[
x[tlll [state[
x[t+l111
return prob
if
‘ main
name
電:
print "
EconomicDyn
担n
i
c
s8
1ページの"
print "
Exercise4
.
2
.
1
5の答え"
print
悶
パス
X，の確率は"， path prob(}
持引数を指定しなければ
持デフォルトのカーネル、初期分布、パスの
非結果をプリント
上のスクリプトを joint utf.py に保存しておいて、以下のように呼び出して使う。
Vol
.1
0
2
8
No.2
持ー*- coding:utf・8 *
-
持 3期ともリセッションの確率
時いずれの期においても、四 MR":mild recession か
神
叩
R":severe recession
2つの内から占つの状態を取ることができるので、全部で 8つのパスがある
非それをすべて列挙して、確率を求め合計する
非
韓
8つのパスのリスト
m
を作る
非
持 joint_utf.py からマルコフ連鎖の遷移行列 pHや初期分布 psi
井関数 path prob(p，
psi，
x) を読み込む
from joint utf import *
e=["MR"，
"SR"l
持各期に可能な状態
m=日
時空のリスト m を作る
for i in e
:
持 for文にはコロン:を忘れない
持
繰り返し部分は半角のスペース
非
4つ分インデン卜すること
for
in e
:
非したがって for丈を入れ子にするには
半角のスペース 4つ分インデント
#
for k in e
:
j，
kl) 井
m.append([
i，
リストクラスのメソッド append
持を使ってリスト [
i，
j
，k
l を m に追加
持
持
韓関数 path_prob(p，
psi，
x
) はパス xの縮率を求める
print
print η 期ともリセッションであるパスとその確率を表示"
枠
s=O.O
非合計を記録する
を作る
ノ、、
F
J
ふ仕
動
上
m
るト
なス
に j
e の
3
ふ品
L
品晶宵
叩ス
、，
ヒ
El
o
可ム﹄
EE
非
ι
0.0 を代人することで
持
for t in m:
B
Pythonを使った経済分析について
持
print
print"パス四， t，刊の確率は"， path prob(x=t)
持 p，
psi にはデフォルト値を使い
x
持
だけを
持パス
s=s+path prob(x=t)
t
t
に変化させる
の確率を s に加える
print
print
printη 期ともリセッションの確率は 8つのパスの確率の和田， s
print
>>>import joint r utf
と実行すると以下のような出力が得られる。
3期ともリセッションであるパスとその確率を表示
パス ['MR'， 'MR'， 'MR'] の確率は 0.1210568
パス['MR'， 'MR'， 'SR'] の確率は 0.0119812
パス ['MR'， 'SR'， 'MR'] の確率は 0.0078232
パス ['MR'， 'SR'， 'SR'] の確率は 0.0075768
パスド SR'， 'MR'， 'MR'] の確率は 0.2371344
パスド SR'， 'MR'， 'SR'] の確率は 0.0234696
パス['SR'， 'SR'， 'MR'] の確率は 0.l499616
パスド SR'， 'SR'， 'SR'] の確率は 0.l452384
3期ともリセッションの確率は
sつのパスの確率の和
joint utf.py を単独で使うにはターミナルの
0.704242
Sプロンプトに対して
$ python joint r utf.py
Economic Dynamics 81 ページの
Exercise 4.2.15 の答え
NG'， '
MR" '
NG')の確率は 0.00084l
パス(，
2
9
.1
0
Vol
30
i
f
NO.2
n
a
m
e =
=
' main 勺ブロックのプリント文が処理されていることが分かる。
S
t
o
c
h
a
s
t
i
cDynamicP
r
o
g
r
a
m
m
i
n
gを解く
初期状態 X E Sからスタートしたときに、政策グを取ったとき得られる満足の合計叫仰
は以下の式で表される。
[
さ
山
一
哨
)
)
]
時):=E
X
(
1
)
が与えられたとき、すべての政策の集合三から、 U
σ(
x
)を最大にする政策 σ*ε 三を見つ
ける問題を考える。この問題を解く方法として
1
.v
a
l
u
ei
t
e
r
a
t
i
o
n
2
.p
o
l
i
c
yi
t
e
r
a
t
i
o
n
の 2つの方法がある。
a
l
u
ef
u
n
c
t
i
o
nVσ仰が以下の (
3
)式で表される
いずれの方法も、次式で定義される v
B
e
l
l
m
a
n方程式を満たすことから派生したものである。
♂(
x
):=sup{v
σ(
x
):σε ~}(x ε S)
(
2
)
これは問題の答の一部である。せめてこれでも分かればと思って、解答をちらっと見ると現
在の状態が xであるとき、 v
a
l
u
ef
u
n
c
t
i
o
nv
*仰は次の B巴
I
lman方程式を満たすと書かれてい
S
q
δ
a
z田町、、、
ε
同
，
，
，〆
11
tj
〆
z
'
a
z
z
、
よい、
+
*
α
qt
ο
r
、
、
，
，
ノ
α
+
Z
U
，、，，，、
xz
r
a
EJ
α
a
u 可l
mG
一
一
本
Z
、
、
，
，
ノ
U
BZ
るO
両辺に関数ザμj が出て来ている、関数を未知とする方程式だ。実は B
e
l
l
m
a
nO
p
e
r
a
t
o
rT を
次のように定義するとき、
(
B
、
Tり (
x
)= m_ax~ U(x-a
)十 p).
v
(α+z
)ゆ(
z
)>(
xε5
)
ぽ r
(x
)l
.~
v
a
l
u
ef
u
n
c
t
i
o
n V可x)は B
e
l
l
m
a
nO
p
e
r
a
t
o
rTの不動点になっており、任意の関数
(
4
)
V(x)に繰り返
I
lmanO
p
e
r
a
t
o
rT を作用させ、とても大きな回数作用させると、 V*仰に限りなく近づ
し B巴
B
a
n
a
c
hの不動点定理)。このことを使って v
a
l
u
ef
u
n
c
t
i
o
n ザ付)
いていくことが知られている (
を求めるのが、 k
u
r
t
z
b
e
l
l
m
a
n
_
r
e
v
i
s
e
d
.
p
yである。この中で T
(
りという関数が定義されている。
o
l
i
c
y も求められるので、
これは関数が U であるときの最大報酬を与える。それと同時に p
e
t
u
mするように変更した。スクリプトの書き方自体はあまり良くないかも
それもついでに r
Pythonを使った経済分析について
知れない。
持 Filename: kurtzbellman rivise.py
非
Author: John Stachurski
非 Date: December 2008
神 Revised:且ugust 2010:Junichi Suga
非
Corresponds to: Listing 5.1
beta，rho，B，M
0.5，0.9，10，5
S
range{B + M + 1
)
井 State space
O，
.
.
.
，B+M
Z
range{B + 1
)
持 Shock space
O，
.
.
.
，B
def U{c):
"
.Utility function・
return c**beta
def phi{z):
.Probability mass function，uniform distribution..
Z
) if 0 <= z <= B else 0
return 1.0 / len(
def Gamma{x):
.The correspondence of feasible actions..
return range{min{x， M) + 1
)
def T{v):
"岡田
An implementation of the Be11man operator.
Parameters: v is a sequence representing a function on S
.
Returns: Tv，a list....
Tv
[
]
PV
[
]
for x in S
:
幹
Compute the value of the objective function for each
非 a in Gamma{x)，and store t
he result in vals
3
1
3
2
Vol
.10
No.2
口
vals
for a in Gamma(x):
y
U(x - a
) + rho * sum(v[a + z]*phi(z) for z in Z
)
vals.append(y)
持 Store the maximum reward for this x in the list Tv
Tv.append(max(vals))
Pv.append(vals.index(max(vals)))
return [Tv，
Pv]
また政策 sigmaを与えたらそのときの報酬を計算する関数 value_oLpoJ
icy関数も retumで返
す v_sigmaの形をよくするために変更を加えた。
神 Filename: kurtzvsigma revise.py
非且uthor: John Stachurski
非 Date: December 2008
持 Revise:August 2010:Junichi Suga
時
Corresponds to: Listing 5.2
from numpy import zeros，dot，array
from kurtzbellman import S， rho， phi， U 非 From listing 5.1
def value_of_policy(sigma):
"Computes the value of following policy sigma."
神 Set up the stochastic kernel p sigma as a 2D array:
N
len(
S
)
p sigma
zeros((N，N))
for x in S
:
for y in S
:
]
p sigma[x，y
phi(y - sigma[x])
非 Create the right Markov operator M sigma:
M sigma
)
lambda h: dot(p sigma，h
Pythonを使った経済分析について
持 Set up the function r sigma as an array:
r sigma
array([U(x - sigma[x]) for x in S])
神 Reshape r sigma into a column vector:
r sigma.reshape((N， 1
)
)
神r sigma
非工 nitialize v sigma to z
ero:
v sigma
zeros(N)
v sigma=array(v sigma)
持 Initialize the discount factor to 1
:
discount
自
l
for i in range(10000):
v sigma
v sigma + discount *r sigma
r sigma
M sigma(r sigma)
discount
discount *rho
return v sigma
そして、
valuei
t
e
r
a
t
i
o
nと p
o
l
i
c
yi
t
e
r
a
t
i
o
nを行うスクリプトを
tesCkurtzbellman_revise.pyとして作成した。
from pylab import plot，
subplot，
show
ones
from numpy import zeros，dot， array，
from kurtzbellman revise import *
from kurtzvsigma rivise import value of po1icy
w=zeros(len(S))
it=O
err=1
while err > 0.000001:
[
v，
p]=T(w)
sa=[
]
3
3
Vol
.10
34
for i i
n range{O，
len{v)・1
):
sa.append{abs(v[i]-w[i]))
err=max{sa)
w=v
subp1ot(2，
2，
1
)
p1ot(w)
2，
3)
subp1ot(2，
p1ot{p)
it=it+l
神show(
)
print 'va1ue function:va1ue iteration'，
it
print{w)
print(p)
持show(
)
韓
q=zeros{len{S))
err=1.0
it=O
whi1e err>O:
print i
t，
q
門司
{
)
E司晶
司ム
C
Y
o
一巾目-
s--1﹂
p
、
，
，w
ヂ
ふ
o， t
可ム司
aM=
U1J
aprL
wrLS
v'=
=va
for i in range(O，
len{pl)-l):
sa.append{abs(pl[i]-q[i]))
err=max(sa)
q=pl
2，
2
)
subp1ot{2，
p10t(
q
)
2，
4
)
subp1ot{2，
plot{v)
it=it+l
teration'
print 'po1icy function:po1icty i
No.2
Pythonを使った経済分析について
3
5
print(q)
w
)
print(
show(
)
以下に実行画面のコピーを示す。かなり不親切なスクリプトで、スクリプトファイルを解読
しないと何をしているか分からない。最初は valueiterationで valuefunction を求めている。
4
0回の繰り返し計算が行われた。 valuefunction と policyfunction を
計算が終了するまで、 1
表示した。次に、 policyiterationで policyfunctionを求めている。最初の policyとして、何に
対しでも Oを選ぶ policyを仮定して計算を初めて 3回の繰り返しで policyfunction が求めら
れている O
va1ue function:va1ue iteration 140
20.749444678896232，
[19.017393871327354， 20.017393871327354， 20.431607433700449，
21.040772645460926， 21.308721837892048， 21.544789815391837， 21.769273465167302，
21.982695230450691， 22.188234882605943， 22.384496450887355， 22.578069018229165，
22.761082924138552， 22.943758737820154， 23.115331613073963， 23.277609273242341]
]
[
0， 0， 0， 0， 1， 1， 1， 2， 2， 3， 3， 4， 5， 5， 5， 5
po1icy function:po1icty iteration
。[
o
.o
.o
.o
.o
.o
.o
.o
.o
.o
.o
.o
.o
. O. o
. 0.1
1 [
0， 0， 0， 0， 1， 1， 1， 1， 2， 2， 2， 3， 3， 4， 4， 4
1
2 [
0， 0， 0， 0， 1，
工
， 1， 2， 2， 3， 3， 4， 4， 5， 5， 5
1
3 [
0， 0， 0， 0， 1， 1， 1， 2， 2， 3， 3， 4， 5， 5， 5， 5
1
繰り返しの様子を図示したものが下図である O 左半分が valueiterationで右半分が policy
i
t
e
r
a
t
i
o
nである。
Vol
.1
0
36
No.2
5
4
3
2
1
2 4
s
5
5
8 10 12 14 16
DO
2 4
6
8 10 12 14 16
6
8 10 12 14 16
24
4
F
fj
41
'-
31
2
1
2 4
6
l
;
;
8 10 12 14 16
2 4
7 Sage:OpenSourceMathematicsSoftwareをf
吏う
Math巴m
a
t
i
c
aや MATLAB を使って演習をしているとき、 3回目が終わった位のときに、
必ず聞かれることは、「卒業しでも Mathematciaや MATLAB を使い続けることができます
か。」ということだが、 Mathematicaや MATLAB を買えば使い続けられるよ、ただし xx万
円だけどねというと、急速に興昧を失い、やる気もうせ、それと同時にヤメる口実を子に入
れることになる。そこで Mathematicaと同じような機能をもつもので wxMaximaというソフ
トだと無料で使える。それに MATLAB と同じ位の機能を持つ Octave もあるからそれも使
えるとアドバイスすると、使ってみると Mathematicaの n
o
t
e
b
o
o
kや MATLABマルチウイン
ドウに慣れているので、物足りないという。世の中お金を出して時間とサービスを買うか自
干の出し具合を適度に調整するものなのであとは、本人次
分で汗を流すかを考えて、お金と i
第ということになる。
P
y
t
h
o
nをターミナルから使ったり、 IDLEで、対話的に使ったり IDLEのエディタで編集し
たりしたけれど、 Winodws7や SnowLeopard，Ubuntunなどで、インストールするだけでは、
日本語を使うところで使い勝手が悪くなる場合も出てくる
6
6WIndows7に標準でついてくる P
o
w
e
r
S
h
e
l
lは、日本語のエンコーデイングがデフォルトでは S
h
i
f
tJ
I
Sなの
でそのままでは使えない。もちろんソフトを使うユーザーインターフェースの所は仕方がないとして、直接
ユーザーが触れない部分では日本語を使わないことで、この問題は解決できるのだと思う
A
Pythonを使った経済分析について
3
7
IPython のような使い勝手で、、いや Mathematicaのような環境で Python を使えないかググ
る(Google) と
、 Sage が検索の網に引っかかる。 Ubuntul0.4上で使えるバイナリがあったの
で
、ためらうことなくインストールした。というのも Ubuntul0.4がすでに Windows上で動
いていたからだ。 7
欝一
品関⋮
i
時
同額一
一鎚一
一時
輯
1
.
.
0
(
'!
o
.
輔 J日.
;
s
a
g
世.
5
.
3
!
j
o
i
n
t時t
f
.
p
x
')
E
.
c
o
n
時昭島良輔l
e
s昌1
ページの
E
<
.
.
c
i
.
.4.2.l5
の寄え
1の確率健在韓鶴4
1
パ A ('総¥'輯¥'総'
l
蝿 d
(
'品窃陣 f
j
畏/措辞4
.
5
.
3
/
j
o
i叫
r立.
t
l
.
o
:
l
l
3態とも早 t"ションで島るパス聖護示
r
[['慣を酎掛金]，
準君安輯" '
S
R
'
]，('鵬官 '
S
R
'，'献]， ('輯
'
S
R
'
.'
S
R
'
J，[
'
S
R
'，'甑
掛'
J，里" 会議有害 '
5
畏
'
J
.[
'3
f
t
'，'
S
R
'
，事
誕R
'
，
]
['SR¥ ，
s
震
¥'主君']J
ψ
r
s
Eeon
時 i
cUyn
喜劇拍車1
ページ母
hereis
e
ぇ2
.
1
5
<
1
.
>
番
え
1<1.>確率は
パス f髭車二金額¥'総'
註酷畢8
4
1
パス亡髄¥'輯二'輔さ]お確率誌邑 1
2
1事部車
パ ^ ['損害櫨ヘ '
S
R
'J<1.>譲率謀長音 1
1
9
8
1
2
パス['輯" 'SR¥'MR']<1.>韓寧は事総'
7
8
2
3
2
パス
f輯九字額¥香 S
R
'] 母護率低華簡157
鵠
パス古事訴¥金援重¥'酷ゴ<1.>確率誌阜 2
3
7
1お 4
安T
7 お金がいらないようにという学生の要望に応えて利用環境を構築するために、 Windows上で簡単にLi
nux
が使えるように仮想環境を整えることにした。まず無料であり USBが使えることで、 V
i
r
t
u
a
l
Bo
xをインス
トールし、さらに比較的使いやすい日本語環境が整った Ubuntuをゲスト OSとして利用していたメンテ
ナンスをしっかりやっていたので、 Ubuntul0.4が動いている環境が既にあった。 Ubuntuの G
u
e
s
t
A
d
d
i
t
i
o
n
s
を再度インストールすると本当にシームレスに H
o
s
tOSである Windowsと G
u
e
s
tOSである Ubuntuとの間
をマウス操作ができて、 Ubuntuを画面いっぱいに大きくして、使っていると、 Windows7の windowと重な
り合って、どちらの OSの windowなのか分からなくなってしまうこともあるくらいだ。また共有フォル夕、
を作っておくと、 Windows7からも Ubuntuからも共有フォルダに読み書きできるのでファイルを共有でき
巴s
tOSの Ubuntuで管理者権限で作ると、平のユーザーは読み書きできなくなるので、
て便利だ。ただし Gu
平のユーザーを所有者としてフォルダを作ることが必要。 H
o
s
tの Windows7側では平のユーザーとしてログ
インしているだろうから、アクセスしやすいデスクトップに共有フォルダを作るとよい。ただ l台の USB
接続のプリンターを 2つの OSで使おうとすると、そのままでは Ubuntuに乗っ取られるようだ。 Windows7
のデスクトップは MacOSXの Dockのようによく使うアプリのアイコンを並べておいて起動できるので便
利である。
3
3
8
Vol
.1
0
NO.2
s
a
g巴を対話モードで使うときには、 p
y
t
h
o
nのファイルをそのまま i
m
p
o
r
t して使えるが、
o
t
e
b
o
o
kで使うときは、 l
m
p
o
r
tコマンドを使うとエラーとなる。した
次の図に示すように n
がって l
o
a
dコマンドを使っている。このあたりのことをもっと調べる必要がある。しかし
n
o
t
e
b
o
o
kを表示するための webs
e
r
v
e
rが起動しており、 n
o
t
e
b
o
o
kのテキスト部分は h
t
m
lで
書けば良いので遠く離れたネット上に置いてリモートで使える等使い方は幅広さを感じる。
8 おわりに
Python を新たに独学の形で習得することを試してみた。他のプログラミング言語を何も
知らない訳ではないので、それを割りヲ I
lミて考えないといけないが、 Pythonは習得が比較
的簡単なのではないだろうか。まず型を宣言しなくても、代人により型を決めてくれる。型
が動的に与えられるので、初心者には簡単である。そして型が決まると、オブジェクト指向
なのでその型のオブジェクトが使えるメソッドが自動的に使えるようになる。これは私たち
がこの世界にあるあり方によく似ている
O
オブジェクトを観察して何ができるのかを知る。
この世にある人工物以外は、世界にすべてビルトインされている。人工物(機械や建物ある
いは制度)ですら、歴史の中ですでに与件として存在しており、その中に投げ込まれる私た
ちのあり方を反映している o Pythonが提供する世界でオブジェクトを探し学習する。経済
分析に使うスクリプトもオブジェクトであり、その可読性の高さがその利用や変更の自由度
を決める。
i
n
u
xのプログラムは、ソースが公開されているので、自分の利用環境
通常私たちが使う L
に対応するバイナリーが見つからなくても、自分でコンパイルしてインストールすることが
できるので、とても便利で、ある。例えば、 P
y
t
h
o
n
2
.
7を試してみたいが、バイナリーがない
ときは、ターミナルからソースをダウンロードして、適当なディレクトリに展開して、イン
a
i
l
s を試したくなっても同じことだ、 R
a
i
l
s
3
.
0のソースファイル
ストールできる o RubyonR
.9
.
2だって自分でソースからインストール
を取ってきて、インストールすると良い。 Ruby1
できる。ただ、 R
a
i
l
s
2と遣うので、インターネットで調べるときには R
a
i
l
s
3のものだけを参
考にしないといけな l'
0また Sageの buildには気をつけた方がよ Lミo sage4.5.3をソースか
u
i
l
d しようとしたとき、 b
u
i
l
dに要する時間の目安として 1時聞から 1
4日と書いている。
らb
よく見ると、 m
o
b
i
l
e端末で b
u
i
l
dすると 1
4日だと書いている。いま使用しているのが、 5
年程前の機械なので、どれくらいかかるのだろうと思ったが、ソースファイルを展開してで
きた s
a
g
e
4
.
5
.
3ディレクトリで makeと入力するだけで簡単だという記述に引きずられて、休
日の朝の 1
0
:
2
0にリターンキーを押してしまった。ターミナルの表示が変わっていくのを、
他の仕事をしながら見ていたが、さすがに 2時間経っとたまらず昼食をとる。その後も他の
仕事をしていたが、 1
6
:
2
0頃になると、モジュールを書き出し、インストールしている。
P
y
t
h
o
nを使った経済分析について
39
b
u
i
l
dが終わったのは、 1
7
:
0
0前だった。 6時間 40分の b
u
i
l
dの後、既存の P
y
t
h
o
nファイルを
動かして動作チェックしたらうまくいった。 Linuxの世界は私のような素人にも使いやすく
なっている
O
分からないことがあれば、ググる (
G
o
o
g
l
e
) ことで解決できる。その際、イン
ターネットから入手できる情報を処理するための c
o
r
e となる知識(その領域の基礎知識と
英語力)の蓄積が大切になると思われる。また学習は自己に働きかけ自己を編集する過程で
あって、この自己を編集する技術やモティベーションを磨き続けることも生涯学習にとって
必要となろう。
参考ウェッブサイトおよび参考文献
P
y
t
h
o
nについての情報を得るためには、 h
t
t
p
:
/
/
w
w
w.
p
y
t
h
o
n
.
o
r
g
!を訪ねるのがよい。そのペ
ージの左にあるリンクの中から DOWNLOAD をクリックして Downloadページに移動する
と
、 Python2.7の i
n
s
t
a
l
l
e
rが並んでいるので、その中かから自分が使っている OS用の
i
n
s
t
a
l
l
e
rをダウンロードしさえすれば簡単にインストールできる。 T
h
i
r
dp
a
r
t
yの s
o
f
t
w
a
r
eの
y
h
t
o
n
3
.
x
.
xより P
y
t
h
o
n
2
.
xの方がよいだろう。 MacOSXには Python
対応状況を考えると、 P
がデフォルトでインストールされているが、より新しいものを使いたいなら、 P
y
t
h
o
n
2
.
7を
インストールするとよい。ただ、 t
h
i
r
dp
a
r
t
yのソフトも P
y
t
h
o
n
2
.
7用をインストールした方
がよ L、。ここで使用した NumPyパッケージは h
t
t
p
:
/
/
s
o
u
r
c
e
f
o
r
g
e
.
n
e
t
!
p
r
o
j
e
c
t
s
/
n
u
m
p
y!
f
i
l
e
s
!を見
.5
.
0
p
y
2
.
7のバイナリーがあるのでインストールするの MATLAB品質の
にいくと、 numpy-1
a
t
p
l
o
t
l
i
bもP
y
t
h
o
n
2
.
7用のものを使う方が良いかも知れな Lしただ、
プロット機能をもっ m
この原稿を書いている時点では、 P
y
t
h
o
n
2
.
7用の b
i
n
a
r
yは h
t
t
p
:
/
/
m
a
t
p
l
o
t
l
i
b
.
s
o
u
r
c
e
f
o
r
g
e
.
n
et
!
の
y
t
h
o
n
2
.
7用のバイナリーがない。そこでソースをダウンロ
ダウンロードページにはな L入 P
ードして .
!
c
o
n
f
i
g
u
r
eから始めてインストールしようとすると、色々なエラーメッセージをは
きだす c 多少の辛抱が必要だ。ターミナル上に吐き出されるログからエラーを読み取り、依
存関係やパージョンの制約をクリアする
A
つの工夫として、 s
e
t
u
p
e
xt
.pyをエディタで編集す
るo b
a
s
e
d
i
rの辞書で、 '
d
a
r
w
i
n
'
:
[
]となっているのを '
d
a
r
w
i
n
'
:
[
'
!
u
s
r
!
X
I
I
R
6
'
]とすると、 f
r
e
巴t
y
p
e
2
とp
n
g
.
hに関するエラーが解決できて、 m
a
t
p
l
o
t
l
i
bをインストールできる。
Sagemathは h
t
t
p
:
/
/
w
w
w
.
s
a
g
e
m
a
t
h
.
o
r
g
/からダウンロードできる。ミラーサイトをみると、
アジアではインドと韓国とロシアしかない。まだまだマイナーなソフトウェアだが、既存の
主だったオープソースを統合し、 Webページから使うというアイデアは面白いと思う。リ
モートの馬力のあるコンビュータで計算して、結果を非力なローカルで文章にまとめるとい
う使い方が可能で、ある。
4
0
Vol
.1
0
NO.2
AppendixL
inuxTips
PATH
PATHの確認
echo $PATH
PATHの追加
先頭に追加
export PATH=/Applications/sage:$PATH
最後に追加
export PATH=$PATH:/Applications/sage
Sage をインストールしたディレクトリである $SAGE_ROOT を PATH に追加する。ここで
$SAGE ROOT は $SAGE ROOT/local の部分であり、私の場合は jApplications/sage であ
る
。
シェルを起動するとき自動的に設定するには、 .bash profile あるいは .bashrc を Vlm
や Emacsで編集して、巴xport文を書き込み保存しておけばよい。
ユーザーのホームディレクトリにいるとして、
vim .bash profile[Enter]
で Vlmエディタが起動する。行のはじめの
だけで、左下に [NewF
i
l
e
] と表示されて何も表
示されないときには、 2つのことが考えられる。 1つ日はファイル名の入力ミス
最初のピ
wF
i
l
e
]の前にあるダブルクオー
に間違いはないか、 [N巴
リオド.や中程のアンダースコア
トの中をチェックしてみる
O
o
2つ目は本当にファイルがない o vim を終了するには
:
q
[
E
n
t
e
r
]半角のコロン:を忘れないように。
b
a
s
h
_
p
r
o
f
i
l
eが表示されているとき、おもむろに
うまく .
1
とタイプすると挿入モードにな
り、文字入力を受け付けるようになる。カーソルを移動させて適当な位置に export文を付
け加えて
[
E
s
c
]でコマンドモードに入り :
w
[
E
n
t
e
けで保存されます。そのまま終了するには :
q
[
E
n
t
e
r
]
でターミナルのプロンプ卜に戻る。くれぐれもコロン:を忘れないように。
さてファイルがあるか無 L功、を調べる方法を考えてみよう。まず自分がどこにいるか知
るO
p
w
d
[
E
n
t
e
r
]で自分のいる場所が分かる。
自分のホームディレクトリにいないことがわかったら、まずは自分のホームディレクトリ
に戻る。実は自分のホームディレクトリに戻るのは簡単で、どこに居ょうと cd[Enter]で一
発で戻ることができる。
T
a
b
]の入力補完機能を使うと、 [
T
a
b
] キーを打った途端
ここで、 vim.bash_p まで、打って [
にp
r
o
f
i
l
e まで入力補完されたら、ファイルがあることがわかる o [
T
a
b
] キーを打ってもファ
P
y
t
h
o
nを使った経済分析について
4
1
イル名が補完されて完成されないようではファイルがあるか怪しい。以下ではファイルやデ
s
(エルヱス)コマンドについて調べてみる。
ィレクトリの情報を表示する I
ファイルやディレクトリの情報を表示
I
s (エルエス)
オプション
-a 隠しファイルを含むすべてのファイルを表示する
-d 指定したディレクトリ自身の情報を表示する
1 ファイルの詳細を表示する
・
-0
ファイルタイプに応じて色分けする
-r エントリ名で降順に並び替えて表示する
-t 最新更新日時を最新のものから順に表示する
-u
-A
I
I
J と同時指定で、最終更新日時ではなく最終アクセス日時を表示する
l
a
J とほぼ同じだが、カレントディレクトリと親ディレクトリは非表示
-F
ディレクトリには/、実行ファイルには*、シンボリックリンクには@を付けて出力する
4
サブディレクトリも含めて表示する
I
sa
[
E
n
t
e
r
]で隠しファイルを含むすべてのファイルを表示するが、表示が止まらずスクロー
ルして最後が表示されて最初の方は見えなくなってしまう。これへの対処法は 3つあって、
1.最後が見えるのだから表示順を逆にすると、元の最初が見える。 rオプションを追加す
る
。
I
sa
r
2 パイフ。を使って I
saの出すデータを処理する。ページャ l
e
s
sをつかう。
I
sa
l
l
e
s
sで隠しファイルも含めた全ファイルを 1画面ずつ表示させるために、 l
saの出
力をパイプ￨を使って l
e
s
sに流し込む。実際に上のコマンドを実行すると、画面左下に
:を表示して入力プロンプトの状態になる。次の画面の表示は{スペース・パー]を押す、
bで前のページに戻る
O
ファイルの有無が確認できてターミナルに戻るには、 q でター
ミナルのプロンプトになる。
4
2
Vol
.1
0
NO.2
3パイフ。を使って I
saの出すデータを処理する。
p
rコマンドでテキストを整形して表示
I
sa
l
p
rt4
[
E
n
t
e
r
]で I
saの 2列 (
2段組み)表示を 4列 (4段組み)表示にする。この
tオプションを付けとかないと、日付や時刻、ヘッダが表示され、改ページ処理
とき -
もされてしまう。ためしに、
4
オプションを付けないで、 1ページあたりの行数を 2
0
行というオプションー 1(エル)ー2
0をつけて実行してみる、つまり I
sa
l
p
r412
0
[
E
n
t
e
r
]
とすると、改ページ処理がされて、ヘッダとして日付と時刻ページ数が表示され、一覧
性に欠ける結果となる。これに tオフションをつけて I
sa
l
p
r412
0t
[
E
n
t
e
r
]で実行す
ると、ヘッダも改ページ処理もされない。一覧性を重視した所望の結果がえられた。
複数のオプションを組み合わせることは先に示したが、たとえば
I
sa
1
t[
E
n
t
e
r
]
とすると、隠しファイルを含むすべてのファイルを、最新更新日時が最新のものから順に
詳細表示することができる
O
またパイフ。を使ってデータを処理することができるので 1
sa
1
t
F
l
g
r
e
p"
*
" [
E
n
t
e
r
l と
すると、隠しファイルを含むすべての実行ファイルを、最新更新日時が最新のものから順に
詳細表示することができる。
l
i
a
sを設定して、コマンド
自分の好みの形のファイルリストの表示方法を見つけたら、 a
ラインで l
sを実行する際には常に自分の好みのフォーマットで表示されるようにするのが
いいだろう。これには、例えば-/.bash~profi1e か γ.bashrc に次のような行を追加すれ
ば
、 d
i
r
[
E
n
t
e
r
]でカレントディレクトリに含まれるディレクトリだけを詳細表示させること
ができる。
a
l
i
a
sd
i
r
=日 sA
1
l
g
r
e
p^
d
"
ここで八dは dで始まる文字列と言う意昧である。早速ターミナルで d
i
rと入力するとディ
レクトリだけの詳細一覧が表示される
O
さらに次の行を加えると、 d
t
[
E
n
t
e
r
] とするだけで、ディレクトリの簡単なリストが 5列に
表示されるようになる O
a
l
i
a
sd
t
=
"
1
sF
l
g
r
e
p/
I
p
rt -5"
P
y
t
h
o
nを使った経済分析について
4
3
長く使っているとどんな a
l
i
a
s を設定したか忘れてしまうので、新しい a
l
i
a
s を設定する前に
l
i
a
s
[
E
n
t
e
r
]で a
l
l
i
a
sの一覧を表示させると良い。
は
、 a
Unixプログラムをソースからインストール
1.圧縮されたソースファイルをダウンロードして展開
ファイルを再帰的にダウンロードするのに都合が良い wgetの最新版をソースからイン
ストールすることを試みる。ターミナルで
open ftp://ftp.gnu.org/pub/gnu/wget/wget-latest.tar.gz[Enterl とすると、 Web
browserが聞いてファイルを保存するかどうか聞いてくるので保存を選ぶとダウンロー
ドが始まる。とても簡単だ。
t
pを使ってファイルを取ることを考えてみる。
それでは、ターミナルから f
t
p
[
E
n
t
e
r
] とするとプロンプトが f
t
p
> に変わる。下のセッションをみると
ターミナルで f
わかるように Name:のフロンフトに士tして anonymous[Enter]としてログインしている。
O
Remotes
y
s
t
巴m t
ypei
sUNIX. であり、いままでやったことが活かされる。 Us
i
n
gb
i
n
a
r
y
mod
巴t
ot
r
a
n
s
f
e
rf
i
l
e
s
.
ノf
イナリモードでファイルをトランスブァするといっている
O
早速今いるディレクトリのファイルとディレクトリの状態を知るために l
s コマンドを
送信する o dで始まる行がディレクトリであるので、最初の該当する行の右端を見ると
gnu となっており、所望のディレクトリが含まれていることが確認できた。そこで目的
のファイルのあるディレクトリへと cdgnu/wget[Enter]で移動する。そして [Tab]による
入力補完を使いながら get wget-latest.tar.gz と入力できたので、ファイ jしがある
ことも分かった。そこで [Enter] でファイル転送が始まり 4秒程で完了。 bye または
c
l
o
s
eで接続を切る o Goodbye.で終了となる。
ftp> open ftp.gnu.org
Connected to ftp.gnu.org.
220 GNU FTP server ready.
Name (ftp.gnu.org:xxxxx): anonymous
230-Due to U.S. Export Regulations， all cryptographic software on this
230-site is subject to the following legal notice:
Vol
.1
0
44
No.2
230230-
This site inc1udes pub1ic1y avai1ab1e encryption source code
230・
which， together with object code resu1ting from the compi1ing of
230-
pub1ic1y avai1ab1e source code，may be exported from the United
230-
States under License Exception "TSU" pursuant to 15 C.F.R. Section
230・
740.13(e).
230230 This 1ega1 notice app1ies to cryptographic software on1y. P1ease see
睡
230-the Bureau of Industry and Security (www.bxa.doc.gov) for more
230・ information about current U.S. regu1ations.
230 Login successful.
Remote system type is UNIX.
Using binary mode to transfer fi1es.
ftp> ls
229 Entering Extended Passive Mode (111433061)
150 Here comes the directory 1isting.
-rw-r--r幽・
20
-rw-r--r--
10
-rw-r--r--
10
-rw-r--r--
10
。
。
。
。
。
。
drwxrwxr-x 302 0
1003
8192 Sep 07 20:35 gnu
drwxrwxr-x
30
1003
4096 Sep 23 2004 gnu+1inux-dis
-rw-r--r--
10
-rw-r--r--
10
drwxr-xr-x
30
lrwxrwxrwx
ユD
drwxr-xr-x
55 0
1rwxrwxrwx
10
drwxr-xr-x
20
drwxr-xr-x
20
drwxr-xr-x
20
1rwxrwxrwx
10
rw-r--r--
工 0
岨
。
。
。
。
。
。
。
。
。
8 Aug 20 2004 CRYPTO.README
17864 Oct 23 2003 MISSING-FILES
4178 Aug 13 2003 MISSING-FILES
1765 Feb 20 2007 README
405121 Oct 23 2003 before-2003・ 08
161919 Sep 19 10:41 find.txt.gz
90 Feb 16 1993 lpf
.README
335164 Sep 19 10:41 ls-lrRt.txt
.gz
4096 Apr 20 2005 mirrors
11 Apr 15 2004 non-gnu ->
4096 Sep 07 17:59 old-gnu
1 Aug 05 2003 pub ->
4096 Nov 08 2007 savannah
4096 Aug 02 2003 third-party
4096 Apr 07 2009 tmp
Pythonを使った経済分析について
drwxr-xr-x
20
0
-rw-r・-r・
開
10
0
4
5
4096 Aug 24 19:42 video
954 Aug 13 2003 we1come.msg
226 Directory send OK.
ftp> cd gnu/wget
250 Directory successfully changed.
ftp> get wget-latest.tar.gz
local: wget-latest.tar.gz remote: wget-1atest.tar.gz
229 Entering Extended Passive Mode (111327851)
150 Opening BINARY mode data connection for wget-latest.tar.gz (246
100¥屯 1*******************************************************************
**********1 2406 KB 524.08 KB/s 00:00 ETA
226 File send OK.
2464747 bytes received in 00:04 (503.12 KB/s)
ftp> bye
221 Goodbye.
それでは、ダウンロードしたファイルをチェックしてみよう。まず
tar tvzf wget_latest.tar.gzI1ess で中身を確認
wget-1
.
12のソースコードであり、 wget田
1
.
12のディレクトリを頂点としてきちんとした
ディレクトリ構造を持っていることが分かる。
tar xvzf wget latest.tar.gz で展開すると、ディレクトリ wgec1
.1
2が作成される。
2
.付属の configureスクリプトで MakeFileを作る
cdwgec1
.1
2として I
sでファイルを見ると、 Makefile.am とか Makefile
i
.
n とかはあるが、
Makefileはない。
s とすると今度は Make ファイルが作ら
そこで、 /configureスクリプトを実行し、再度 I
れている。
3
.makeでコンノ tイル
4
.makei
n
s
t
a
l
lでインスト -)1>
Vol
.1
0
4
6
これで私の利用環境
No.2
CPUの種類や OSなどに適合するようにカスタマイズした wgetの
最新版をインストールすることができた。
wgetを使う
wgetをつかってカレントディレクトリからリンクを張っているファイ lしをすべてダウン
ロードする。 MATLABや P
y
t
h
o
n、R などのスクリプトを探していると、ファイルのリンク
がずらっと表示されるページに遭遇することがある。一つ一つ別名でリンク先を保存と地道
にやれば良いが、時間とやる気が無いときには wgetを使うと良い。
例えばスクリプトのファイルへのリンクが並んだページの URL(UniformR
e
s
o
u
r
c
eL
o
c
a
t
o
r
)
がh
t
t
p
:
/
/
w
w
w.
x
x
x
x
.
x
xx/yyyy/
c
o
d
e
/であるとすると、
ターミナルで wget rn
p ・11 ・且 pyhttp://www.xxxx.xxx/yyyy/code/ とすると、
帰的にリンクを辿ることを指定し、 n
p で上位のディレクトリ yyyy
パラメーター r で再J
へのリンクは辿らないことを指定する。念のためにリンクを辿る深さを -11 として l回に
制限する。間違っても o としてはいけない。無限にリンクを追いかけてしまうので。さら
に
、
-A
pyとして、ファイルの拡張子が pyであるファイルのみに限定する。すると 40個程
のファイルがあっという聞にダウンロードされる O
wget[OPTION]…[URL]
-B
[ベース URLJ
-i ファイル名]
ベース URLを指定する。
指定されたファイルの URLを順に処理する
-B
司
o [ファイル名]
-T
[タイムアウト時間]
O
でベース URLを指定すれば相対 URLでよい。
指定したファイル名で取得結果を保存する
タイムアウト時聞を秒で指定する。

Download Report