pythonのclassを使って、1農家の生産行動と1消費者の消費行動をざっくりシミュレーションしてみる②　～結果を見る回～

前回、経済モデルを使って条件をどのように設定したか説明しましたので、今回はどんなコードでその状況を構築したのか説明します。
条件設定から振り返りたい方は、前回の[pythonのclassを使って……①](https://qiita.com/kuru_kirino/items/cb3bb8e4176804fdbf28)の方をご確認ください。

##今回 やることの確認
①前回設定したモデルのパラメータに具体的な数値を代入する。
②コードを確認する。
③結果を確認する。

##①パラメータに具体的な数値を代入する
###他の経済シミュレーションにおけるパラメータの扱い
先日、読んだ経済セミナー増刊号（2016）の「進化する経済学の実証分析」によると、
最近のマクロ経済学系の論文では、経済モデルが現実を表せているかどうかを計量経済的に検証した後、推定したパラメータをその経済モデルに代入してから反実仮想シミュレーション(もしも補助金が交付されなかったら農家の生産はどうなっていたか？　とか)を行う事がしばしばあるようです。
前回の記事で僕が紹介したTaylor(2005)とDyer(2006)もメキシコの農村調査から生産関数や効用関数（コブ＝ダグラス型）の中のパラメータを簡単に推定してシミュレーションしていました。
一方、このシミュレーションではパラメータの推定は行わず、適当に選んだ数値を代入しました。
以下、今回のシミュレーションでパラメータにあてはめた数値の紹介をします。

###農家の生産関数
```math
\Pi = p_A*Q_A+p_B*Q_B+200*L_n...(1) \\ Q_A = 3*L_A^{0.3} *K_A^{0.7}...(2) \\ Q_B = L_B^{0.7} *K_B^{0.3}...(3) \\ 24 = L_A+L_B+L_n...(4) \\ 10 = K_A+K_B...(5)
```
$p_A(Aの初期価格)$:100
$p_B(Bの初期価格)$:300
※上記価格は時間に沿って変化する。

###消費者の効用関数
```math
U = c_A^{0.4}+c_B^{0.6}...(1) \\ I_t = I_{t-1} + 5000 - C_{t-1} > C_t = p_A*c_A + p_B*c_B...(2:予算制約) \\ Q_{A,t} + s_{A,t-1} > c_{A,t}...(3:供給+余り＞消費) \\ Q_{B,t} + s_{B,t-1} > c_{B,t}...(4:供給+余り＞消費) \\ c_A + c_B < 40...(5:満腹制約)
```
$p_A(Aの初期価格)$:100
$p_B(Bの初期価格)$:300

###その他（価格形成と次期の生産物持ち越し）
①次期の予算：$I_{t+1} =I_t + 5000 - c_t$
②次期の持ち越し量：$s_{t+1} = s_t + Q_t - c_t$
③次期の価格決定方法
　if $(Q_t + s_{t-1} - c_t)/(Q_t + s_{t-1}) < 0.3$:
　    $p_{t+1} = 1.01 * p_t$
　else:
　    $p_{t+1} = 0.99 * p_t$

なんでこのように設定したのかは前の記事で説明しました。
そちらをご確認ください。

##②コード
```python
"""
sell&buy.pyプログラム
農家と消費者が生産と購入を繰り返す
農家は利潤を最大化し、消費者は効用を最大化する
農家は自家労働力と自家資本を生産に用いる
次期の価格は、前期におけるストック＋生産量－消費量を基に変動する
"""

#　モジュールのインポート
import random
import numpy as np
import matplotlib.pyplot as plt

#　グローバル変数
alpha = 0.3   #　生産物Aの生産にかかわるパラメータ
beta = 0.7    #　生産物Bの生産にかかわるパラメータ
gamma = 0.4  #　消費者が各生産物を消費するときに得る効用にかかわるパラメータ
w = 200       #　賃金
appe = 40     #　1人の人間が食べられる最大量を40とした
SEED = 32768  #　シード値
R = 100        #　取引の繰り返し回数

#　クラス定義 Farmer
class Farmer:
    """農家を表現するクラスの定義"""
    def __init__(self, cat):   #　コンストラクタ
        self.category = cat
        self.maxprofit = 0   #　利潤の最大値（ひとまず0）
        self.mqa = 0          #　利潤最大時のAの生産量
        self.mqb = 0          #　利潤最大時のBの生産量
        self.mLA = 0         #　以下、利潤最大時の各生産要素の使用量
        self.mKA = 0
        self.mLB = 0
        self.mKB = 0
        self.mLn = 0
        self.L = 24           #　自家労働量
        self.K = 10           #　自家資本量
    def solvePmax(self):   #　利潤を最良化するようにAとBを生産する
        """利潤の最良値をリセットする"""
        self.maxprofit = 0
        """利潤の最適化問題を解く"""
        for i in range(nlimit):
            self.calcprofit(Market)   #　利潤が計算される（profitが計算される）
    def calcprofit(self, Market):   #　利潤計算
        """生産要素の自家制約を満たす組み合わせをランダムに得る"""
        LA = int(random.randrange(24))       #　０から２４以内の整数を選ぶ
        LB = int(random.randrange(24 - LA))  #　０から２４－LA 以内の整数を選ぶ
        Ln = self.L - LA - LB                #　余った労働を農外労働に投入
        KA = int(random.randrange(10))
        KB = self.K - KA
        """生産関数"""
        qa = 3 * LA ** (alpha) * KA ** (1 - alpha)
        qb = 1 * LB ** (beta) * KB ** (1 - beta)
        """利潤関数"""
        profit = Market.pA * qa + Market.pB * qb + w * Ln
        """最良解を更新"""
        if profit >= self.maxprofit:   
            self.maxprofit = profit
            self.mqa = qa
            self.mqb = qb
            self.mLA = LA
            self.mKA = KA
            self.mLB = LB
            self.mKB = KB
            self.mLn = Ln
#　クラスFarmerの定義終わり

#　クラス定義 Consumer
class Consumer:
    """消費者を表現するクラスの定義"""
    def __init__(self, cat):   #　コンストラクタ
        self.category = cat
        self.maxutility = 0   #　効用の最大値（ひとまず0）
        self.mca = 0     #　生産物Aの消費量
        self.mcb = 0     #　生産物Bの消費量
        self.C = 0       #　総支出額
        self.I = 5000    #　初期状態の予算として貯蓄S=5000を与えた
    def solveUmax(self):   #　効用を最良化するようにAとBを消費する
        """効用の最良値と予算制約をリセットする"""
        self.maxutility = 0
        I = self.I + 5000   #　一回当たりの収入は5000
        """効用の最適化問題を解く"""
        for i in range(nlimit):
            self.calcutil(I, Market)   #　効用が計算される（utilityが計算される）
    def calcutil(self, I, Market):   #　効用計算
        """ランダムに生産物AとBの消費量を選ぶ"""
        ca = int(random.randrange(appe))       #　０からappe=40以内の整数を選ぶ
        cb = int(random.randrange(appe - ca))  #　０からappe=40 - ca 以内の整数を選ぶ
        """今期の支出総額を計算"""
        C = Market.pA * ca + Market.pB * cb
        """効用関数"""
        utility = ca ** (gamma) + cb ** (1 - gamma)  
        """最良解を更新"""
        if I > C and Market.sa > ca and Market.sb > cb:
            if utility >= self.maxutility:   #　最良解を更新
                self.maxutility = utility
                self.mca = ca
                self.mcb = cb
                self.C = C
                self.I = I - C   #　次期の予算を計算するので、支出を引いた
            
#　クラスConsumerの終わり        

#　クラス定義 Market
class Market:
    def __init__(self):
        self.pA = 100   #　初期時点の農産物価格
        self.pB = 300   
        self.sa = 2 * (12 ** 0.3) * (5 ** 0.7)   #　自家保有する生産要素の半分をそれぞれAとBの生産に振り分けた場合の値を用いた
        self.sb = 1 * (12 ** 0.7) * (5 ** 0.3)
    def save(self, f):   #　農家が生産した分が貯蔵庫に入れられる
        for i in range(len(f)):   
            self.sa += f[i].mqa
            self.sb += f[i].mqb  
    def eaten(self, c):   #　消費者が消費した分が貯蔵庫から惹かれて、更に1/3だけネズミに食われる
        for i in range(len(c)):   #　人間に食わて次期の貯蔵が減る
            self.sa -= c[i].mca
            self.sb -= c[i].mcb
        self.sa *= 2 / 3          #　ネズミに食われて次期の貯蔵が減る
        self.sb *= 2 / 3
    def price(self, c):   # 次期の価格を計算
        Tca = 0
        Tcb = 0
        for i in range(len(c)):
            Tca += c[i].mca
            Tcb += c[i].mcb
        if (self.sa * 3 / 2) / (self.sa * 3 / 2 + Tca) <= 0.3:
            self.pA *= 1.01
        else:
            self.pA *= 0.99
        if (self.sb * 3 / 2) / (self.sb * 3 / 2 + Tcb) <= 0.3:
            self.pB *= 1.01
        else:
            self.pB *= 0.99               
        
#　下請け関数の定義
#　次期の計算を統括するcalcn()関数の定義
def calcn(t, f, c, Market):
    #　①時間、前期の価格、前期の生産物余りをリストに入れた
    pAlist.append(Market.pA)   #　価格
    pBlist.append(Market.pB)
    salist.append(Market.sa)   #　生産物の余り
    sblist.append(Market.sb)
    tlist.append(t + 1)        #　時間
    #　②農家の利潤最良化
    for i in range(len(f)): 
        f[i].solvePmax()
    #　③市場に出回る生産物の生産量やその他の要素の集計値
    mqa = 0
    mqb = 0
    maxprofit = 0
    mLA = 0
    mKA = 0
    mLB = 0
    mKB = 0
    mLn = 0
    for i in range(len(f)):   
        mqa += f[i].mqa
        mqb += f[i].mqb
        maxprofit += f[i].maxprofit
        mLA += f[i].mLA
        mKA += f[i].mKA
        mLB += f[i].mLB
        mKB += f[i].mKB
        mLn += f[i].mLn    
    #　④各農家アウトプットの集計値をリストに入れる
    maxprofitlistF.append(maxprofit)   
    mqalistF.append(mqa)
    mqblistF.append(mqb)
    mLAlistF.append(mLA)
    mKAlistF.append(mKA)
    mLBlistF.append(mLB)
    mKBlistF.append(mKB)
    mLnlistF.append(mLn)
    #　⑤市場に農産物が出回る
    Market.save(f)
    #　⑥消費者の効用最良化
    for i in range(len(c)):   
        c[i].solveUmax()
    #　⑦市場で起こる生産物の消費量やその他の要素の集計値  
    mca = 0
    mcb = 0
    maxutility = 0
    C = 0
    mI = 0
    for i in range(len(c)):   
        mca += c[i].mca
        mcb += c[i].mcb
        maxutility += c[i].maxutility
        C += c[i].C
        mI += c[i].I    
    #　⑧各消費者アウトプットの集計値をリストに入れた
    maxutilitylistC.append(maxutility)   
    mcalistC.append(mca)
    mcblistC.append(mcb)
    ClistC.append(C)
    IlistC.append(mI)
    #　⑨次期に持ち越される生産物の余り
    Market.eaten(c)
    #　⑩次期の価格
    Market.price(c)      
    

#　メイン実行部
#　初期化
random.seed(SEED)      #　乱数の初期化
Market = Market()      #　変数にクラスMarketのインスタンスを入れる
#　リストの中にカテゴリ0のFarmerとConsumerのインスタンスを複製生成（する予定。今回は一人ずつ）。なおリストでインスタンスを作ると、変数単体では使用できないようだ。
f = [Farmer(0)]
c = [Consumer(0)]
#　農家の状態のリスト化
maxprofitlistF = []
mqalistF = []
mqblistF = []
mLAlistF = []
mKAlistF = []
mLBlistF = []
mKBlistF = []
mLnlistF = []
#　消費者の状態のリスト化
maxutilitylistC = []
mcalistC = []
mcblistC = []
ClistC = []
IlistC = []
#　価格と生産物の余りのリスト化
pAlist = []
pBlist = []
salist = []
sblist = []
#　時間のリスト
tlist = []
#　試行回数の入力
nlimit = int(input("最良解を追及する回数を指定してください。大きい数を指定する程、最適解に近づきます。："))

#　エージェントシミュレーション
for t in range(R):
    calcn(t, f, c, Market)   #　次時刻の状態を計算
#　シミュレーションパートの終わり

#　特定のリストにおける数値の小数点以下調整
#　利潤、生産量、効用、C、s、収入を調整する
maxprofitlistF = [round(maxprofitlistF[n]) for n in range(len(maxprofitlistF)) ]
mqalistF = [round(mqalistF[n], 1) for n in range(len(mqalistF))]
mqblistF = [round(mqblistF[n], 1) for n in range(len(mqblistF))]
maxutilitylistC = [round(maxutilitylistC[n], 1) for n in range(len(maxutilitylistC))]
ClistC = [round(ClistC[n]) for n in range(len(ClistC))]
IlistC = [round(IlistC[n]) for n in range(len(IlistC))]
salist = [round(salist[n], 2) for n in range(len(salist))]
sblist = [round(sblist[n], 2) for n in range(len(sblist))]


#　グラフの表示
print("利潤","\n",maxprofitlistF,"\n","\n","A生産量","\n",mqalistF,"\n","\n","B生産量","\n",mqblistF,"\n","\n","LA","\n",mLAlistF,"\n","\n","LB","\n",mLBlistF,"\n","\n","Ln","\n",mLnlistF,"\n","\n","効用","\n",maxutilitylistC, "\n","\n","A消費","\n",mcalistC,"\n","\n","B消費","\n",mcblistC,"\n","\n","C","\n",ClistC,"\n","\n","pA","\n",pAlist,"\n","\n","pB","\n",pBlist,"\n","\n","sa","\n",salist,"\n","\n","sb","\n",sblist,"\n","\n","収入","\n",IlistC)
pA = plt.plot(tlist, pAlist)
pB = plt.plot(tlist, pBlist)
plt.title("the price change of product A and B")
plt.xlabel("time")
plt.ylabel("price")
plt.legend((pA[0], pB[0]), ("pA", "pB"), loc=4)

plt.show()
#　sell&buy.py
```

###全体の構造
####classを使った
題名にも書いてある通り、classを使ってみました。
ここで作ったクラスは3種類あり、「Farmer」、「Consumer」、「Market」です。
Farmerは利潤を最大化するようにAとBを生産し、
Consumerは制約下で効用を最大化するようにAとBを消費し、
MarketはAとBの在庫管理と価格調整を行います。

####最適化問題の解き方
他の経済シミュレーションにおいて、最適化問題がどのように解かれているのか分かりませんが、
ここでは指定回数だけランダムに選択変数の組み合わせを試し、効用か利潤が一番大きくて制約条件を満たしているときの組み合わせをどんどん更新していく形をとりました。

####コードを書く上で参考にした本
コードの書き方は、小高知宏（2018）：『Pythonによる数値計算とシミュレーション』に則っていますが、小高さんが使った言葉を誤った意味合いで解釈した部分があるかもしれません。

####エージェントの数が1人なのに、リスト変数にFarmerとConsumerのインスタンスを入れた
この理由は、後々エージェントの数を増やそうと思ったからです。
今後、マルチエージェントっぽくしていく際は、パラメータが違うエージェントを増やそうと思ってます。
エージェントのパラメータをばらけさせる際に、
①確率分布に従ってパラメータが微調整される感じにするか、
②パラメータが他の外生的な変数の影響を受けて変化する構造にするか、悩みます。

##③結果の確認
僕がランしてみたところ、結果は以下のように出力されました。
なお、長いので一部省略。

###結果
最良解を追及する回数を指定してください。大きい数を指定する程、最適解に近づきます。：100 
 A生産量 
 [11.7, 21.2, 22.0, 16.3, 11.4, 21.2, 17.2, 17.9, 19.4, 14.4, 17.9, 17.2, 23.9, 19.4, 15.8, 22.6, 19.5, 19.4, 15.8, 15.8, 19.4, 19.5, 21.2, 17.2, 17.9, 22.0, 14.6, 22.6, 16.3, 17.2, 19.5, 19.4, 25.0, 17.8, 19.5, 23.9, 23.9, 25.0, 17.9, 27.0, 22.6, 14.4, 19.5, 21.2, 19.5, 19.4, 23.9, 15.8, 17.2, 26.1, 17.9, 19.4, 21.2, 17.2, 21.2, 26.1, 19.4, 22.6, 26.1, 22.6, 23.9, 21.2, 23.9, 22.6, 19.0, 26.1, 25.0, 27.9, 19.5, 22.6, 21.2, 25.0, 19.5, 22.6, 22.0, 21.2, 19.4, 27.0, 25.0, 27.0, 24.0, 19.4, 27.9, 21.2, 22.6, 19.5, 27.9, 27.0, 25.0, 27.0, 20.8, 30.1, 25.0, 23.9, 23.9, 26.1, 29.4, 19.4, 17.2, 25.0] 
 
 B生産量 
 [5.4, 0.0, 3.8, 4.9, 6.3, 1.6, 1.0, 3.2, 0.0, 2.3, 3.2, 2.2, 1.0, 2.6, 3.2, 1.6, 2.0, 1.0, 3.2, 2.0, 2.6, 2.7, 2.2, 0.0, 2.0, 1.2, 5.9, 1.6, 4.3, 1.0, 3.2, 1.0, 1.6, 3.0, 1.2, 1.6, 2.6, 2.6, 1.2, 1.6, 2.2, 3.0, 5.7, 2.6, 1.2, 3.1, 1.6, 3.2, 1.6, 0.0, 1.2, 1.6, 2.6, 0.0, 2.6, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 1.6, 1.0, 2.3, 2.6, 0.0, 2.2, 3.8, 2.2, 1.0, 2.2, 1.2, 1.0, 3.8, 1.0, 3.9, 0.0, 3.1, 1.0, 2.0, 2.2, 3.1, 1.0, 1.6, 3.2, 1.0, 1.0, 1.6, 0.0, 1.2, 1.0, 0.0, 2.6, 2.2, 1.6, 0.0, 2.2, 1.0, 1.6] 
 
 LA 
 [1, 4, 6, 3, 2, 4, 2, 3, 3, 2, 3, 2, 6, 3, 2, 5, 4, 3, 2, 2, 3, 4, 4, 2, 3, 6, 3, 5, 3, 2, 4, 3, 7, 4, 4, 6, 6, 7, 3, 9, 5, 2, 4, 4, 4, 3, 6, 2, 2, 8, 3, 3, 4, 2, 4, 8, 3, 5, 8, 5, 6, 4, 6, 5, 5, 8, 7, 10, 4, 5, 4, 7, 4, 5, 6, 4, 3, 9, 7, 9, 8, 3, 10, 4, 5, 4, 10, 9, 7, 9, 5, 13, 7, 6, 6, 8, 12, 3, 2, 7] 
 
 LB 
 [7, 0, 5, 6, 7, 2, 1, 4, 0, 2, 4, 3, 1, 4, 4, 2, 2, 1, 4, 2, 4, 3, 3, 0, 2, 1, 7, 2, 5, 1, 4, 1, 2, 3, 1, 2, 4, 4, 1, 2, 3, 3, 9, 4, 1, 5, 2, 4, 2, 0, 1, 2, 4, 0, 4, 0, 1, 1, 0, 0, 0, 1, 2, 1, 2, 4, 0, 3, 5, 3, 1, 3, 1, 1, 5, 1, 7, 0, 5, 1, 2, 3, 5, 1, 2, 4, 1, 1, 2, 0, 1, 1, 0, 4, 3, 2, 0, 3, 1, 2] 
 
 Ln 
 [16, 20, 13, 15, 15, 18, 21, 17, 21, 20, 17, 19, 17, 17, 18, 17, 18, 20, 18, 20, 17, 17, 17, 22, 19, 17, 14, 17, 16, 21, 16, 20, 15, 17, 19, 16, 14, 13, 20, 13, 16, 19, 11, 16, 19, 16, 16, 18, 20, 16, 20, 19, 16, 22, 16, 16, 20, 18, 16, 19, 18, 19, 16, 18, 17, 12, 17, 11, 15, 16, 19, 14, 19, 18, 13, 19, 14, 15, 12, 14, 14, 18, 9, 19, 17, 16, 13, 14, 15, 15, 18, 10, 17, 14, 15, 14, 12, 18, 21, 15] 
 
 A消費 
 [23, 11, 28, 16, 12, 14, 19, 17, 21, 7, 23, 16, 23, 20, 16, 12, 26, 17, 16, 7, 20, 21, 19, 19, 15, 21, 10, 22, 10, 23, 16, 22, 20, 20, 16, 26, 20, 24, 16, 25, 23, 9, 13, 28, 19, 20, 20, 6, 22, 24, 20, 19, 20, 20, 12, 28, 14, 25, 26, 1, 39, 21, 13, 29, 18, 24, 20, 32, 17, 24, 20, 26, 15, 17, 25, 23, 11, 23, 20, 34, 16, 21, 22, 18, 27, 11, 33, 23, 28, 28, 28, 28, 28, 15, 25, 27, 27, 10, 22, 26] 
 
 B消費 
 [12, 1, 3, 5, 5, 2, 1, 3, 0, 2, 3, 2, 1, 1, 3, 2, 2, 0, 2, 3, 1, 1, 4, 0, 0, 2, 6, 1, 4, 1, 2, 2, 1, 1, 1, 2, 3, 2, 1, 2, 1, 3, 6, 2, 1, 3, 1, 2, 2, 0, 0, 2, 3, 3, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 2, 3, 0, 2, 3, 2, 1, 2, 1, 1, 4, 1, 2, 1, 3, 1, 2, 1, 3, 1, 2, 3, 0, 1, 2, 2, 2, 2, 2, 0, 2, 1, 0, 2, 0, 2] 
 
 C 
 [5900, 1414, 3700, 3101, 2724, 2048, 2244, 2660, 2185, 1341, 3310, 2299, 2753, 2459, 2668, 1930, 3438, 1859, 2391, 1707, 2520, 2651, 3363, 2161, 1723, 3042, 3006, 2861, 2394, 2977, 2486, 3209, 2696, 2717, 2250, 3795, 3368, 3608, 2330, 3797, 3276, 2066, 3479, 4192, 2764, 3521, 2946, 1405, 3502, 3200, 2693, 3165, 3627, 3627, 1681, 3884, 2255, 3468, 3642, 141, 5463, 3256, 2145, 4387, 3135, 4317, 2915, 5279, 3389, 4184, 3326, 4563, 2614, 2947, 5039, 3890, 2334, 3889, 4052, 5613, 3134, 3588, 4385, 3176, 4869, 2676, 5205, 3972, 5115, 5115, 5115, 5115, 5115, 2537, 4867, 4959, 4704, 2351, 3833, 5166] 

 収入 
 [4100, 7686, 8986, 10886, 13162, 16114, 18870, 21210, 24025, 27684, 29374, 32075, 34322, 36863, 39195, 42265, 43827, 46968, 49578, 52870, 55350, 57700, 59337, 62176, 65452, 67410, 69404, 71543, 74149, 76171, 78685, 80477, 82781, 85064, 87813, 89018, 90650, 92042, 94712, 95916, 97639, 100573, 102095, 102902, 105138, 106617, 108672, 112267, 113764, 115565, 117872, 119707, 121080, 121080, 124399, 125515, 128260, 129792, 131150, 136009, 135546, 137290, 140145, 140758, 142623, 143307, 145392, 145112, 146723, 147539, 149212, 149649, 152036, 154089, 154050, 155160, 157827, 158937, 159885, 159272, 161138, 162550, 163164, 164989, 165120, 167444, 167239, 168267, 168152, 168152, 168152, 168152, 168152, 170616, 170748, 170789, 171085, 173734, 174901, 174735]

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/611511/34089639-49ed-505d-edf4-65c6c5246c79.png)


###結果を見た感想
本気で考察結果をまとめようとするともっと長くなるので、結果を見た感想をざっくりと述べたいと思います。

####①見せ方（図表）
matplotlibを使った結果の出力が結構難しいことがわかりました。
今回、価格の変化だけ図にして、他の変数の変化は適当にリストを貼っただけで見せましたが、全ての変数を棒グラフにしてもよかったなぁ。次の課題です。

####②価格の変化
なぜ価格の変化だけ表にしたかというと、価格が均衡価格みたいなところに収束するのか気になったからです。
図を見てみると、
生産物A（価格が低くて資本集約的な穀物をイメージした）の価格がどんどん上がったことが分かります。そして、Aの消費量はほとんど一定で、生産量はじわじわと増加したように見えます。
価格が上がったAを農家がより多く作るのは普通ですが、消費者がAを消費する量が変化していない点は不思議です。
おそらく理由は、満腹度制約があるせいで、収入の増加量に対して消費額の増加量が過少になってしまい、予算が増えすぎて予算制約が意味をなさなくなったからかもしれません。
その場合、消費の配分は、予算と関係なく効用関数のパラメータによって定まる消費の均衡状態（価格が変化しないときにおける、一番好きな消費の配分）に収束するかもしれないです。
うーん、奥深い、、、。
満腹制約は外した方がいいのかな。もしくは、エージェントを増やしたらよくなるのかな、、、

####③生産要素の投入傾向
生産要素の投入傾向を一言で説明すると
「農外労働の賃金が高すぎて、みんな副業目的でしか農業をやらない」です。
Lnにほとんどの時間が投資されたことがわかりましたが、賃金200という設定は大きすぎたのかもしれないです。
また、生産物B（価格が高くて労働集約的な果物をイメージ）がほとんど生産されず、Aだけがたくさん作られました。この条件設定だと、農家は少し初期価格が高いBではなくて、生産要素の投入生産性が高いAを多く作る傾向があったということが言えます。なんか、日本の兼業米農家みたいですね。
ちなみに、パラメータに違う数値を入れると全く違う結果が出力されます。やっぱりパラメータの値にもう少しこだわる必要があるのかもしれないです。

####⑤世界は複雑なんだなぁと思った
友達が言っていたのですが、セルオートマトンの設計者は、最後にどのような模様が描かれるのか、プログラムをランするまで分からないそうです。
今回、僕（このプログラムの神）は農家と消費者の行動を模倣できそうなモデルを作ってみたわけですが、ランする前はどんな結果が出るか分からなかったし、結果を見た後になんでその結果になったのかの理由すら断定できませんでした。つくづく、科学者は無謀な挑戦をしているのかもしれないと感じます。現実はこのシミュレーションよりもっと複雑なわけですからね。
結局、人間にできることは、「①何らかの現象を説明する際に骨となる基礎理論を見つける。②特定の現象を説明する際にはいろいろな骨をくっつけまくって肉にする」ことだけなんだと思います。


ということで、唐突ですがここで終わります。
次は自給的農家のシミュレーションを作るか、もしくはpygameを勉強して農家がワラワラ動く感じの簡単なゲームを作ろうかなと思っています。


##参考文献
～論文～
[1]Alain de Janvry, Marcel Fafchamps and Elisabeth Sadoulet(1991) “Peasant household behavior with missing markets: some paradoxes explained” Economic Journal. Vol.101, pp.1400-1417.
[2]George A. Dyer, Steve Boucher, and J. Edward Taylor(2006) “Subsistence response to market shocks” American journal of agricultural economics. Vol. 88, pp. 279-291.
[3]J.Edward Taylor, George A. Dyer, Antonio Yu'nez-Naude(2005) "Disaggregated Rural Economywide Models for Policy Analysis" World Development. vol. 33, pp. 1671-1688

～テキストブック～
[1]小高知宏（2018）：『Pythonによる数値計算とシミュレーション』, オーム社
[2]Jake VanderPlas(2016) "Python Data Science Handbook: Essential Tools for Working with Data", O'Reilly Media