pythonのclassを使って、自給農家の生産・消費行動を超ざっくり簡単にシミュレーションしてみる② コードと結果を見る

じゃあ、かったるいおさらいは省いてさっそくコードと結果を見てみましょう。
条件設定の説明は①をご確認ください。

##コード

"""
subsistance.pyプログラム
自給的農家が生産をする
自給的農家はランダムに生産計画を繰り返し練り、効用が最大化される計画を実行する
自給的農家は自家労働力と自家資本を生産に用いる
1~10期：初期価格で効用最大化を行う
11~20期:Bの価格が1度上昇して、その価格のまま効用最大化を行う
21~30期:oの価格が1度低下して、その価格のまま効用最大化を行う
31~40期:Bの価格が1度上昇し、尚且つoが1度低下して、その価格のまま効用最大化を行う
"""

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

#　グローバル変数
alpha = 0.3   #　生産物Aの生産にかかわるパラメータ
beta = 0.7    #　生産物Bの生産にかかわるパラメータ
gamma1 = 0.2   #　消費者が自給食料を消費するときに得る効用にかかわるパラメータ
gamma2 = 0.45  #　消費者がその他の財（テレビとかNETFLIXとか）を消費するときに得る効用にかかわるパラメータ
w = 200       #　賃金
SEED = 1234  #　シード値
R = 40        #　繰り返し回数の設定

#　クラス定義 
class Autarky:
    """自給的農家を表現するクラスの定義"""
    def __init__(self, cat):   #　コンストラクタ
        self.category = cat
        self.utility = 0      #　最大効用（ひとまず0）
        self.co = 0           #　効用最大時のOの消費量
        self.l = 0            #　効用最大時の余暇時間
        self.C = 0            #　効用最大時の総消費額
        self.I = 5000         #　効用最大時の予算(5000としておく)
        self.profit = 0       #　効用最大時の利潤
        self.qa = 0           #　効用最大時のAの生産量
        self.qb = 0           #　効用最大時のBの生産量
        self.LA = 0           #　以下、効用最大時の各生産要素の使用量
        self.KA = 0
        self.LB = 0
        self.KB = 0
        self.Ln = 0
    def solveUmax(self):   #　利潤を最良化するようにAとBを生産する
        """効用の最良値をリセットする"""
        self.utility = 0
        """利潤をランダムに計算した後、その利潤に基づく効用の最適化問題を解く"""
        for i in range(nlimit):
            self.calcprofit(Market)                   #　利潤が計算される
            I = self.I + self.profit                  #　農家の収入は農業労働と農外労働による収入の和
            self.calcutility(I, Market)               #　効用が計算される
    def calcprofit(self, Market):    #　利潤計算
        """生産要素の自家制約を満たす組み合わせをランダムに得る"""
        LA = int(random.randrange(18))                #　０から18以内の整数を選ぶ
        LB = int(random.randrange(18 - LA))           #　０から18－LA 以内の整数を選ぶ
        Ln = int(random.randrange(18 - LA - LB))      #余った労働を農外労働に投入
        l = 18 - LA - LB - Ln
        KA = int(random.randrange(10))
        KB = 10 - KA
        """生産関数"""
        qa = 3 * LA ** (alpha) * KA ** (1 - alpha)
        qb = 1 * LB ** (beta) * KB ** (1 - beta)
        """利潤関数"""
        profit = Market.pB * qb + w * Ln
        """情報を更新""" 
        self.profit = profit
        self.qa = qa
        self.qb = qb
        self.LA = LA
        self.KA = KA
        self.LB = LB
        self.KB = KB
        self.Ln = Ln
        self.l = l
    def calcutility(self, I, Market):   #　効用の最大値を計算する
        """ランダムにoの消費量を選ぶ"""
        maxco = I / Market.po
        co = int(maxco)
        """今期の支出総額を計算"""
        C = Market.po * co
        """効用関数"""
        utility = self.qa ** (gamma1) + co ** (gamma2) + self.l ** (1 - gamma1 - gamma2)
        """効用の最良解を更新"""
        if utility >= self.utility:   #　最良解を更新
            self.utility = utility
            self.co = co
            self.C = C
            self.I = I - C   #　次期の予算を計算するので、支出を引いた        
#　クラスAutarkyの定義終わり     

#　クラス定義 Market
class Market:
    def __init__(self):
        self.pB = 100   #　初期時点の農産物価格
        self.po = 300   
    def price_B(self):   #　pBが50円値上げ
        self.pB += 50
    def price_o(self):   #　poが50円値下げ
        self.po -= 50
#　クラスMarketの定義終わり

        
#　下請け関数の定義
#　価格無変化時における次期の計算を統括するcalcn_n()関数の定義
def calcn_n(t, a, Market):
    #　①今期初めにおける価格をリストに入れる
    pBlist.append(Market.pB)   #　価格
    polist.append(Market.po)
    tlist.append(t)            #　時間
    #　②自給的農家の効用最良化
    for i in range(len(a)):
        a[i].solveUmax()


#　pBが1度変化した時における次期の計算を統括するcalcn_B()関数の定義
def calcn_B(t, a, Market):
    #　①pB価格の変化
    Market.price_B()
    #　②今期初めにおける価格をリストに入れる
    pBlist.append(Market.pB)   #　価格
    polist.append(Market.po)
    tlist.append(t)            #　時間
    #　③農家の効用最良化
    for i in range(len(a)):
        a[i].solveUmax()

    
#　poが1度変化した時における次期の計算を統括するcalcn_o()関数の定義
def calcn_o(t, a, Market):
    #　①po価格の変化
    Market.price_o()
    #　②今期初めにおける価格をリストに入れる
    pBlist.append(Market.pB)   #　価格
    polist.append(Market.po)
    tlist.append(t)            #　時間
    #　③農家の効用最良化
    for i in range(len(a)):
        a[i].solveUmax()


#　pBとpoが1度変化した時における次期の計算を統括するcalcn_Bo()関数の定義
def calcn_Bo(t, a, Market):
    #　①pBとpo価格の変化
    Market.price_B()
    Market.price_o()
    #　②今期初めにおける価格をリストに入れる
    pBlist.append(Market.pB)   #　価格
    polist.append(Market.po)
    tlist.append(t)            #　時間
    #　③農家の効用最良化
    for i in range(len(a)):
        a[i].solveUmax()


#　要素の集計値計算及びリストへの追加を行うCalc_list()関数の定義
def calc_list(a):
    qa = 0
    qb = 0
    profit = 0
    LA = 0
    KA = 0
    LB = 0
    KB = 0
    Ln = 0
    l = 0
    cb = 0
    co = 0
    utility = 0
    C = 0
    I = 0
    for i in range(len(a)):
        qa += a[i].qa
        qb += a[i].qb
        profit += a[i].profit
        LA += a[i].LA
        KA += a[i].KA
        LB += a[i].LB
        KB += a[i].KB
        l += a[i].l
        Ln += a[i].Ln
        co += a[i].co
        utility += a[i].utility
        C += a[i].C
        I += a[i].I
    utilitylistA.append(utility)
    colistA.append(co)
    ClistA.append(C)
    IlistA.append(I)
    profitlistA.append(profit)
    qalistA.append(qa)
    qblistA.append(qb)
    LAlistA.append(LA)
    KAlistA.append(KA)
    LBlistA.append(LB)
    KBlistA.append(KB)
    LnlistA.append(Ln)
    llistA.append(l)


#　メイン実行部
#　初期化
random.seed(SEED)      #　乱数の初期化
Market = Market()      #　変数にクラスMarketのインスタンスを入れる
#　リストの中にカテゴリ0のAutarkyのインスタンスを複製生成（する予定。今回は一人ずつ）
a = [Autarky(0)]
#　農家の状態のリスト化
utilitylistA = []
colistA = []
ClistA = []
IlistA = []
profitlistA = []
qalistA = []
qblistA = []
LAlistA = []
KAlistA = []
LBlistA = []
KBlistA = []
LnlistA = []
llistA = []
#　価格と生産物の余りのリスト化
pBlist = []
polist = []
#　時間のリスト
tlist = []
#　試行回数の入力
nlimit = int(input("最良解を追及する回数を指定してください。大きい数を指定する程、最適解に近づきます。："))

#　シミュレーション
for t in range(R):
    if t <= 10:                     #　1~10期：初期価格
        calcn_n(t, a, Market)
        calc_list(a)
    if t == 11:                     #　11期：B価格変化
        calcn_B(t, a, Market)
        calc_list(a)
    if (t > 11) and (t <= 20):      #　12~20期：B価格変化後の安定
        calcn_n(t, a, Market)
        calc_list(a)
    if t == 21:                     #　21期：o価格変化
        calcn_o(t, a, Market)
        calc_list(a)
    if (t > 21) and (t <= 30):      #　22~30期：o価格変化後の安定
        calcn_n(t, a, Market)
        calc_list(a)
    if t == 31:                     #　31期：Bとo価格変化
        calcn_Bo(t, a, Market)
        calc_list(a)
    if (t > 31) and (t <= 40):     #　32~40期：Bとo価格変化後の安定
        calcn_n(t, a, Market)
        calc_list(a)
#　シミュレーションパートの終わり

#　重要な変数のリストをそのまま表示
print("利潤",profitlistA,"\n","\n","A生産",qalistA,"\n","\n","B生産",qblistA,"\n","\n","LA",LAlistA,"\n","\n","LB",LBlistA,
      "\n","\n","Ln",LnlistA,"\n","\n","l", llistA, "\n","\n",, "効用",utilitylistA, "\n","\n", "o消費",colistA,"\n","\n","C",ClistA,
      "\n","\n","pA",pBlist,"\n","\n","po",polist,"\n","\n", "収入",IlistA)

#　sell&buy.py

##結果
今回は、最良解を追求する回数を10000にしました。シード値は1234で固定です。

###価格

設定した通り、10期・20期・30期でBとoの価格が増減しました。

###AとBの生産量

Aの変化がグダグダw
10・20・30期を境にU字でAの生産量が変化しました。なんでだ、、、？w
Bの生産量は割と一定なように見えますね。

###A,o,lの消費量

oの消費量がじわじわと増加してきています。価格が減少傾向にあるため、当然でしょう。
一方、l（余暇）の消費量は一定ですね。

###生産要素の投入量

みにくいw。積み上げ棒グラフにすればよかったですね。まあいいか。
何となくですが、時系列的な負の自己相関傾向があるような気がします。前期の投入量が大きいとき、次期は下がってますよね。なんだこれ。

###利潤と消費額

まあ、一定ですね。
換金作物価格は増加傾向にあるので、利潤は増加するのかと思いましたが、自給穀物の生産増加が邪魔して利潤が上がってないですね。単純に解釈するなら、自給生産しているせいで農家の収入が低位に安定してしまったということかな。これが自給生産の罠なのかな。

###効用

おーーww
やっぱり、消費できる量自体は増えてるわけだから、効用は増加してますね。
特に20から上がってるから、このシミュレーションでは、換金作物Bの価格増加よりも、嗜好品oの価格減少の方が自給農家的には良いイベントだったってことですね。
なるほどー。

##問題点
###最大化問題
前回指摘しましたが、最大化問題のやり方がよくないです。というのも、10000回ランダムに生産要素の投入配分を変えまくって、一番いい結果を採用する形で最良解を見つけているため、シード値を変えると結果がまるで変ったりします。利潤最大化を目的にする場合は問題ないのですが、利潤の算出をまたいで効用最大化しようとすると安定しなくなります。
やっぱり、数式的に最大化問題を解く形にしないといけないのかな。紙と鉛筆で解く方法はわかるのですが、どういうコードを使えばいいのかな。
また、なんとなくイメージですが、最尤法によって回帰式のパラメータを推定するときにニュートン＝ラフソン法が使われますが、あーいう感じでじわじわと解に迫っていくアプローチができないかなとも考えてます。
このままでもいい気もしますけどね。

###関数系はコブ＝ダグラス型で良いのか？　また、パラメータは一定でいいのか？
何となくですが、農家は農法に習熟するにつれて生産性を上げていく気がします。この場合、パラメータも改善されるべきなのではないでしょうか？　もしくは、全要素生産性が上がるべきなのではないでしょうか？　この点が気になります。
また、関数型をコブ＝ダグラス型に固定していますが、規模の経済性やパラメータの変動も考慮するために、もっと複雑な関数系にすべきではないでしょうか？　さて、どうしようか。

###そろそろマルチエージェントにしない？
一応、これでクラスは3種類用意できるようになりました。「農家」「消費者」「自給的農家」です。
また、趣味的な家庭菜園農家もしくは自分の農法にやりがいを感じている有機農家などもモデル化するアイデアがあります。これらの人たちをぶっこんだ一つの世界を作りたいなーと思ってます。

以上です。
次は、趣味的な家庭菜園農家or自分の農法にやりがいを感じてる有機農家の生産モデルをシミュレーションなしでここに挙げてみようかなと思ってます。
そのあと、マルチエージェントシミュレーションを走らせてみるかな。
最大化問題のやり方は、ひとまずこのままでいきます。