発展途上国における「農家の最適生産（ミクロ）」と「市場の価格形成」を連関させたシミュレーション＿②コード

1, 導入

じゃあ、コードお見せしていきます。
内部の理論に関する説明は、前記事「～～①」を参照でお願いします。

そういえば、これを言っておかないとダメな条件があったのを忘れていたので、冒頭で言っておきます。
①農家は4人いる
②一人は絶対にnet-sellerになるが、残り3人は最適な形態（seller or buyer）を選ぶ。
③保有農地面積にグラデーションがある（20ha, 20ha, 20ha, 20ha）
です（ここ、あとでコードを見返して気づいたんですが、農地面積が揃ってました。ということは、全員同じ行動をするはずでは、、、？）。

2, 環境

python 3.8.3
scipy 1.5.0
statsmodels 0.12.0
tqdm 4.50.2
をjupyter notebook上で使ってます。

3, コード

※randomというモジュールをimportし忘れていたので、動かなかったです。修正しました（このバージョンでは、randomをimportしなくても良かったですけど、、）。

#　モジュールのインポート
import numpy as np
import pandas as pd
from scipy.optimize import minimize
from tqdm import tqdm_notebook as tqdm
import statsmodels.api as sm
import random

#　グローバル変数（4か月を1サイクルとする）
ao = 0.06287930254251227                      #　その他の財を消費するときの係数(0.06)
al = 0.7171948682330835                      #　余暇を消費するときの係数(0.72)
aoo = -0.01                       #　その他の財消費における2乗の係数
aLl = -0.01                       #　余暇消費における2乗の係数
aol = 0.15                       #　まんべんなく使えるときのうれしさ
w = 75.25114451327885                         #　賃金(時給)
tw = 20                          #　労働市場で労働を売る際にかかる単位当たり取引費用
ta = 5024207380 / 8324248200 / 10 #　Aを売るときの取引費用（初期価格の10分の1と簡単に設定）
r = 2000 * 4                     #　地価（4か月一括で借りる前提の地代）
SEED = 12345                     #　シード値
R = 20                           #　繰り返し回数の設定
L = (24 - 6) * 30 * 4            #　初期賦存労働力（4か月あたり時間）
Yield = 5024207380 / 396397.92   #　単収（kg / ha）
Equip = 1 / 52                   #　土地装備率（ha / 時間）
appe = 0.25 * 30 * 4             #　満腹制約（4か月当たりに食べる量）
pAlist = [5024207380 / 8324248200] #　作物価格の初期値（peso/kg）
milk_meat = 130/4 + 25/4         #　食費

#　クラス定義 
class Farmer:
    """自給的農家を表現するクラスの定義"""
    def __init__(self, cat):   #　コンストラクタ
        self.category = cat
        if cat == "a": #　絶対にnet-sellerになる中規模農家
            self.oriA = 20
        if cat == 0:   #　小規模農家
            self.oriA = 2         #　保有土地面積(ha)
        if cat == 1:   #　中規模農家
            self.oriA = 20         #　保有土地面積(ha)
        if cat == 2:   #　大規模農家
            self.oriA = 200         #　保有土地面積(ha)
        self.oriL = L         #　保有労働時間
        self.utility = 0      #　最大効用（ひとまず0）
        self.co = 0           #　効用最大時のOの消費量
        self.ca = 0           #　効用最大時の穀物購入量
        self.l = 0            #　効用最大時の余暇時間
        self.C = 0            #　効用最大時の総消費額
        self.profit = 0       #　効用最大時の利潤
        self.qa = 0           #　効用最大時のAの生産量
        self.supply = 0       #　効用最大時の市場へのAの供給量
        self.LA = 0           #　以下、効用最大時の各生産要素の使用量
        self.AA = 0
        self.Ln_re = 0        #　受け取る労働力
        self.Ln_gi = 0        #　上げる労働力
        self.An_re = 0
        self.An_gi = 0
        self.expected_pAlist = [] #　Aの期待価格のリスト
        self.type = "net-seller"  #　農家のタイプ
        self.save = 70000         #　初期前期収入(自給農家の収入が70000だったので、大体このくらいにしておく)
    """
    効用最大化問題
    x = [co, La, Ln, l, Aa, An, co_n]
    net-seller→net-buyerの効用最大化問題を解いていき、効用が大きかったら更新する。
    """    
    def solveUmax(self):   #　利潤を最良化するようにAとBを生産する
        if t <= 0:
            pA = pAlist[t]
        else:
            pA = self.expected_pAlist[t-1]
        po = 30
        oriA = self.oriA
        oriL = self.oriL
        income_pre = self.save

        if self.category == "a":   #　絶対にnet-sellerになる人(cat = "a")
            """Seller"""
            """効用関数"""
            def utility(x):   #　二次式モデル、今期の割引率1、来季の割引率0.5で来期の余暇は勘定に入れない
                return -1 * (ao*x[0] + al*x[3] + aoo*x[0]*x[0] + aLl*x[3]*x[3] + aol*x[0]*x[3]) -0.5 * (ao*x[6] + aoo*x[6]*x[6])
            """流動性制約"""
            def cons_liquidity(x):
                return income_pre - po * x[0] - (w + tw) * x[2] - r * x[5] - milk_meat
            """来期の流動性制約"""
            def cons_liquidity_next(x):   #　来期も今期と同じ生産計画になる前提で設定(補助金は1ha当たり1300peso in 2016)
                return (pA - ta) * (Yield * (x[4] + x[5]) - appe) + x[4] * 1300 - po * x[6] - (w + tw) * x[2] - r * x[5] - milk_meat
            """空腹制約"""
            def appetite(x):
                return Yield * (x[4] + x[5]) - appe
            """土地面積と投入労働力の関係性"""
            def labor_land(x):
                return Equip * (x[1] + x[2]) - x[4] + x[5]
            """賦存要素制約"""
            def labor_cons(x):
                return oriL - x[1] - x[3]
            def land_cons(x):
                return oriA - x[4]
            """各変数の符号制約"""
            def co(x):
                return x[0]
            def La(x):
                return x[1]
            def Ln(x):
                return x[2]
            def l(x):
                return x[3]
            def Aa(x):
                return x[4]
            def An(x):
                return x[5]
            def co_n(x):
                return x[6]
            def Anmax(x):
                return 20 - x[5]

            cons = ({"type": "ineq", "fun": cons_liquidity}, {"type": "ineq", "fun": cons_liquidity_next}, {"type": "ineq", "fun": appetite}, 
                    {"type": "ineq", "fun": labor_land}, {"type": "ineq", "fun": La}, {"type": "ineq", "fun": Ln}, {"type": "ineq", "fun": l},
                    {"type": "ineq", "fun": co}, {"type": "ineq", "fun": Aa}, {"type": "ineq", "fun": An}, {"type": "ineq", "fun": co_n}, 
                    {"type": "eq", "fun": labor_cons}, {"type": "ineq", "fun": land_cons}, {"type": "ineq", "fun": Anmax})

            #　Aの生産に要素を全振りしている状態を初期値とする。
            x0 = [0, oriL, 0, 0, oriA, 0, 0]

            res = minimize(utility, x0, constraints=cons, method="SLSQP")   #　最適化問題を解くコマンド

            #　結果が負になったときの対応
            while (res.x[0] < 0) or (res.x[1] < 0) or (res.x[2] < 0) or (res.x[3] < 0) or (res.x[4] < 0) or (res.x[5] < 0) or (res.x[6] < 0):
                x0 = []
                for X in res.x:
                    if X < 0:
                        X = 0
                    x0.append(X)
                res = minimize(utility, x0, constraints=cons, method="SLSQP")

            #　農地投入がかなり大きいのに、労働投入が０になったときの対応(初期値をちょっと変えてやりなおし)
            roop = 0.1
            while (res.x[1] + res.x[2] <= 0) and (res.x[4] + res.x[5] > 1):
                x0 = [0, oriL - roop, 0, 1, Equip * (oriL - roop), 0]
                res = minimize(utility, x0, constraints=cons, method="SLSQP")
                #　結果が負になったときの対応2
                while (res.x[0] < 0) or (res.x[1] < 0) or (res.x[2] < 0) or (res.x[3] < 0) or (res.x[4] < 0) or (res.x[5] < 0) or (res.x[6] < 0):
                    x0 = []
                    for X in res.x:
                        if X < 0:
                            X = 0
                        x0.append(X)
                    res = minimize(utility, x0, constraints=cons, method="SLSQP")
                roop = roop + 0.1

            """計算した結果をしまっていく"""
            self.LA = res.x[1]
            self.Ln_re = res.x[2]
            self.AA = res.x[4]
            self.An_re = res.x[5]
            self.co = res.x[0]
            self.l = res.x[3]
            self.utility = (ao*res.x[0] + al*res.x[3] + aoo*res.x[0]*res.x[0] + aLl*res.x[3]*res.x[3] + aol*res.x[0]*res.x[3])
            self.profit = (pA - ta) * (Yield * (res.x[4] + res.x[5]) - appe)
            self.save = income_pre - po * res.x[0] - (w + tw) * res.x[2] - r * res.x[5] - milk_meat
            self.qa = Yield * (res.x[4] + res.x[5])
            self.supply = Yield * Equip * (res.x[4] + res.x[5]) - appe
            self.C = po * res.x[0] + (w + tw) * res.x[2] + r * res.x[5] + milk_meat


        else:   #　必ずしもnet-sellerにならない人（cat = 数字)
            """Seller"""
            """効用関数"""
            def utility(x):   #　二次式モデル、今期の割引率1、来季の割引率0.5で来期の余暇は勘定に入れない
                return -1 * (ao*x[0] + al*x[3] + aoo*x[0]*x[0] + aLl*x[3]*x[3] + aol*x[0]*x[3]) -0.5 * (ao*x[6] + aoo*x[6]*x[6])
            """流動性制約"""
            def cons_liquidity(x):
                return income_pre - po * x[0] - (w + tw) * x[2] - r * x[5] - milk_meat
            """来期の流動性制約"""
            def cons_liquidity_next(x):   #　来期も今期と同じ生産計画になる前提で設定(補助金は1ha当たり1300peso in 2016)
                return (pA - ta) * (Yield * (x[4] + x[5]) - appe) + x[4] * 1300 - po * x[6] - (w + tw) * x[2] - r * x[5] - milk_meat
            """空腹制約"""
            def appetite(x):
                return Yield * (x[4] + x[5]) - appe
            """土地面積と投入労働力の関係性"""
            def labor_land(x):
                return Equip * (x[1] + x[2]) - x[4] + x[5]
            """賦存要素制約"""
            def labor_cons(x):
                return oriL - x[1] - x[3]
            def land_cons(x):
                return oriA - x[4]
            """各変数の符号制約"""
            def co(x):
                return x[0]
            def La(x):
                return x[1]
            def Ln(x):
                return x[2]
            def l(x):
                return x[3]
            def Aa(x):
                return x[4]
            def An(x):
                return x[5]
            def co_n(x):
                return x[6]
            def Anmax(x):
                return 20 - x[5]

            cons = ({"type": "ineq", "fun": cons_liquidity}, {"type": "ineq", "fun": cons_liquidity_next}, {"type": "ineq", "fun": appetite}, 
                    {"type": "ineq", "fun": labor_land}, {"type": "ineq", "fun": La}, {"type": "ineq", "fun": Ln}, {"type": "ineq", "fun": l},
                    {"type": "ineq", "fun": co}, {"type": "ineq", "fun": Aa}, {"type": "ineq", "fun": An}, {"type": "ineq", "fun": co_n}, 
                    {"type": "eq", "fun": labor_cons}, {"type": "ineq", "fun": land_cons}, {"type": "ineq", "fun": Anmax})

            #　Aの生産に要素を全振りしている状態を初期値とする。
            x0 = [0, oriL, 0, 0, oriA, 0, 0]

            res = minimize(utility, x0, constraints=cons, method="SLSQP")   #　最適化問題を解くコマンド

            #　結果が負になったときの対応
            while (res.x[0] < 0) or (res.x[1] < 0) or (res.x[2] < 0) or (res.x[3] < 0) or (res.x[4] < 0) or (res.x[5] < 0) or (res.x[6] < 0):
                x0 = []
                for X in res.x:
                    if X < 0:
                        X = 0
                    x0.append(X)
                res = minimize(utility, x0, constraints=cons, method="SLSQP")

            #　農地投入がかなり大きいのに、労働投入が０になったときの対応(初期値をちょっと変えてやりなおし)
            roop = 0.1
            while (res.x[1] + res.x[2] <= 0) and (res.x[4] + res.x[5] > 1):
                x0 = [0, oriL - roop, 0, 1, Equip * (oriL - roop), 0]
                res = minimize(utility, x0, constraints=cons, method="SLSQP")
                #　結果が負になったときの対応2
                while (res.x[0] < 0) or (res.x[1] < 0) or (res.x[2] < 0) or (res.x[3] < 0) or (res.x[4] < 0) or (res.x[5] < 0) or (res.x[6] < 0):
                    x0 = []
                    for X in res.x:
                        if X < 0:
                            X = 0
                        x0.append(X)
                    res = minimize(utility, x0, constraints=cons, method="SLSQP")
                roop = roop + 0.1

            """計算した結果をしまっていく"""
            self.LA = res.x[1]
            self.Ln_re = res.x[2]
            self.AA = res.x[4]
            self.An_re = res.x[5]
            self.co = res.x[0]
            self.l = res.x[3]
            self.utility = (ao*res.x[0] + al*res.x[3] + aoo*res.x[0]*res.x[0] + aLl*res.x[3]*res.x[3] + aol*res.x[0]*res.x[3])
            self.profit = (pA - ta) * (Yield * (res.x[4] + res.x[5]) - appe)
            self.save = income_pre - po * res.x[0] - (w + tw) * res.x[2] - r * res.x[5] - milk_meat
            self.qa = Yield * (res.x[4] + res.x[5])
            self.supply = Yield * (res.x[4] + res.x[5]) - appe
            self.C = po * res.x[0] + (w + tw) * res.x[2] + r * res.x[5] + milk_meat


            """Buyer"""
            """x = [co, La, Ln, l, Aa, An, ca, co_n]"""
            """効用関数"""
            def utility(x):   #　二次式モデル、今期の割引率1、来季の割引率0.5で来期の余暇は勘定に入れない
                return -1 * (ao*x[0] + al*x[3] + aoo*x[0]*x[0] + aLl*x[3]*x[3] + aol*x[0]*x[3]) -0.5 * (ao*x[7] + aoo*x[7]*x[7])
            """流動性制約"""
            def cons_liquidity(x):
                return income_pre - po * x[0] - milk_meat
            """来期の流動性制約"""
            def cons_liquidity_next(x):
                return (w - tw) * x[2] + r * x[5] + x[4] * 1300 - po * x[7] - pA * x[6] - milk_meat
            """空腹制約"""
            def appetite(x):
                return Yield * x[4] + x[6] - appe
            """土地面積と投入労働力の関係性"""
            def labor_land(x):
                return Equip * x[1] - x[4]
            """賦存要素制約"""
            def labor_cons(x):
                return oriL - x[1] - x[2] - x[3]
            def land_cons(x):
                return oriA - x[4] - x[5]
            """各変数の符号制約"""
            def co(x):
                return x[0]
            def La(x):
                return x[1]
            def Ln(x):
                return x[2]
            def l(x):
                return x[3]
            def Aa(x):
                return x[4]
            def An(x):
                return x[5]
            def ca(x):
                return x[6]
            def co_n(x):
                return x[7]

            cons = ({"type": "ineq", "fun": cons_liquidity}, {"type": "ineq", "fun": cons_liquidity_next}, {"type": "ineq", "fun": appetite}, {"type": "ineq", "fun": labor_land},
                   {"type": "ineq", "fun": co}, {"type": "ineq", "fun": La}, {"type": "ineq", "fun": Ln}, {"type": "ineq", "fun": l},
                   {"type": "ineq", "fun": Aa}, {"type": "ineq", "fun": An}, {"type": "ineq", "fun": ca},  {"type": "ineq", "fun": co_n}, {"type": "eq", "fun": labor_cons}, 
                   {"type": "ineq", "fun": land_cons})

            #　Aの生産に要素を全部りしている状態を初期値とする。
            x0 = [0, oriL, 0, 0, oriA, 0, 0, 0]

            res = minimize(utility, x0, constraints=cons, method="SLSQP")   #　最適化問題を解くコマンド

            #　結果が負になったときの対応
            while (res.x[0] < 0) or (res.x[1] < 0) or (res.x[2] < 0) or (res.x[3] < 0) or (res.x[4] < 0) or (res.x[5] < 0) or (res.x[6] < 0) or (res.x[7] < 0):
                x0 = []
                for X in res.x:
                    if X < 0:
                        X = 0
                    x0.append(X)
                res = minimize(utility, x0, constraints=cons, method="SLSQP")

            #　農地投入がかなり大きいのに、労働投入が０になったときの対応(初期値をちょっと変えてやりなおし)
            roop = 0.1
            while (res.x[1] <= 0) and (res.x[4] > 1):
                x0 = [0, oriL - roop, 0, 1, Equip * (oriL - roop), 0, 0, 0]
                res = minimize(utility, x0, constraints=cons, method="SLSQP")
                #　結果が負になったときの対応2
                while (res.x[0] < 0) or (res.x[1] < 0) or (res.x[2] < 0) or (res.x[3] < 0) or (res.x[4] < 0) or (res.x[5] < 0) or (res.x[6] < 0) or (res.x[7] < 0):
                    x0 = []
                    for X in res.x:
                        if X < 0:
                            X = 0
                        x0.append(X)
                    res = minimize(utility, x0, constraints=cons, method="SLSQP")
                roop = roop + 0.1

            """効用が大きくなったなら、計算した結果をしまっていく"""
            if - 1 * res.fun > self.utility:
                self.Ln_gi = res.x[2]
                if res.x[4] < 0.0001:   #　投入農地面積があまりにも小さいときは、その人はnet-buyer判定
                    self.LA = 0
                    self.AA = 0
                    self.An_gi = res.x[5] + res.x[4]
                    self.ca = res.x[6] + Yield * res.x[4]
                    self.qa = 0
                    self.profit = (w - tw) * res.x[2] + r * res.x[5] + res.x[4] * 1300
                    self.save = income_pre + (w - tw) * res.x[2] + r * res.x[5] + res.x[4] * 1300 - po * res.x[0] - pA * res.x[6]
                    self.C = po * res.x[0] + pA * res.x[6]
                else:
                    self.LA = res.x[1]
                    self.AA = res.x[4]
                    self.An_gi = res.x[5]
                    self.ca = res.x[6]
                    self.qa = Yield * res.x[4]
                    self.profit = (w - tw) * res.x[2] + r * res.x[5]
                    self.save = income_pre + (w - tw) * res.x[2] + r * res.x[5] - po * res.x[0] - pA * res.x[6]
                    self.C = po * res.x[0] + pA * res.x[6]
                self.utility = (ao*res.x[0] + al*res.x[3] + aoo*res.x[0]*res.x[0] + aLl*res.x[3]*res.x[3] + aol*res.x[0]*res.x[3])
                self.co = res.x[0]
                self.l = res.x[3]
                self.supply = 0
                self.type = "net-buyer"
    #def solveUmax終わり

    def expect_price(self):
        if (1 <= t) and (t <=2):
            self.expected_pAlist.append(pAlist[t] * (1 + (pAlist[t] - pAlist[0]) / pAlist[0]))
        if 3 <= t:
            pAlist_pre = pAlist[:-1]
            pAlist_cur = pAlist[1:]
            df = pd.DataFrame({"pA":pAlist_cur, "pA_pre":pAlist_pre})
            y = df["pA"]
            X = df["pA_pre"]
            ols_res = sm.OLS(y, sm.add_constant(X)).fit()   #　OLSで回帰
            expected_pA = ols_res.predict()[-1]   #　時期のpA予測
            if expected_pA < 0:   #　予測値が負になったら、価格が0になると予想させる
                self.expected_pAlist.append(0)
            else:
                self.expected_pAlist.append(expected_pA)
    #def expect_price終わり

    def Sell(self):
        if self.type == "net-seller":   #　net-sellerだけ農産物を売ったことによる収入がある
            pA = pAlist[-1]             #　均衡価格
            if (pA - ta) * (Yield * (self.AA+ self.An_re) - appe) + self.AA * 1300 > 0:           #　利益が正じゃなかったら売らない(廃棄？)
                profit = (pA - ta) * (Yield * (self.AA+ self.An_re) - appe) + self.AA * 1300
            else:
                profit = 0
            #　時期に使う貯蓄の更新
            self.save = self.save + profit


#　クラス定義 Demand
class Demand:
    def __init__(self):
        self.demand_A = 70000
    def upgrade_demand(self, Supply):   #　次の均衡価格が決まった後における需要量（=その時の供給量）
        self.demand_A = Supply



#　下請け関数の定義
#　農家が将来価格を予想して生産する
def FarmerDecision(a):
    for i in range(len(a) - 1):
        a[i].expect_price()   #　時期価格を予測
        a[i].solveUmax()      #　今期の生産を決行

#　需要と供給を合わせる関数の定義
def Match(a, Demand):
    #　①fixされた生産量の総和を計算
    Supply = 0
    for i in range(len(a) - 1):
        Supply = Supply + a[i].supply
    #　②Supply = Demandとなる価格とその時の需要を決定, 式：Supply = 70000 + 64166 * pA_pre - 64166 * pA_n
    #　基準価格：前期の価格との差で計算してしまうと、前期の価格が高くなればなるほど、Supplyが大きい場合に時期の価格が高くても問題なくなる。
    #　したがって、基準価格は5024207380 / 8324248200(2018年飼料トウモロコシ価格)とした。
    #　だから、式：Supply = 70000 + 64166 * 5024207380 / 8324248200 - 64166 * pA_n
    pA_pre = 5024207380 / 8324248200   #　メキシコの実際のトウモロコシの価格
    pA_n = (70000 + 64166 * pA_pre - Supply) / 64166
    #　③需要の更新と時間・価格のリスト追加
    Demand.upgrade_demand(Supply)
    if pA_n <= 0:
        pAlist.append(0)
    else:
        pAlist.append(pA_n)
    #　④均衡価格に応じた収入をnet-sellerが得る
    for i in range(len(a) - 1):
        a[i].Sell()

#　要素の集計値計算及びリストへの追加を行うCalc_list()関数の定義
def calc_list(a):
    qa = 0
    profit = 0
    LA = 0
    Ln_re = 0
    Ln_gi = 0
    AA = 0
    An_re = 0
    An_gi = 0
    l = 0
    co = 0
    utility = 0
    C = 0
    seller = 0
    buyer = 0
    for i in range(len(a)):
        qa += a[i].qa
        profit += a[i].profit
        LA += a[i].LA
        AA += a[i].AA
        Ln_re += a[i].Ln_re
        Ln_gi += a[i].Ln_gi
        An_re += a[i].An_re
        An_gi += a[i].An_gi
        l += a[i].l
        co += a[i].co
        utility += a[i].utility
        C += a[i].C
        if a[i].type == "net-seller":
            seller += 1
        if a[i].type == "net-buyer":
            buyer += 1
    utilitylistA.append(utility)
    colistA.append(co)
    ClistA.append(C)
    profitlistA.append(profit)
    qalistA.append(qa)
    LAlistA.append(LA)
    AAlistA.append(AA)
    An_relistA.append(An_re)
    An_gilistA.append(An_gi)
    Ln_relistA.append(Ln_re)
    Ln_gilistA.append(Ln_gi)
    llistA.append(l)
    SellerNum.append(seller)
    BuyerNum.append(buyer)


#　メイン実行部
#　初期化
random.seed(SEED)      #　乱数の初期化
Demand = Demand()      #　変数にクラスMarketのインスタンスを入れる
#　リストの中にカテゴリ0のAutarkyのインスタンスを複製生成（する予定。今回は一人ずつ）
a = [Farmer("a"), Farmer(0), Farmer(0), Farmer(0)]
#　農家の状態のリスト化
utilitylistA = []
ca_plistA = []
colistA = []
ClistA = []
profitlistA = []
qalistA = []
qblistA = []
LAlistA = []
AAlistA = []
An_relistA = []
An_gilistA = []
Ln_relistA = []
Ln_gilistA = []
llistA = []
SellerNum = []
BuyerNum = []
#　価格と生産物の余りのリスト化
polist = []
#　時間のリスト
tlist = []

#　シミュレーション
for t in tqdm(range(R)):
    if t <= 10:
        tlist.append(t)     #　tをtlistに追加
        FarmerDecision(a)   #　農家が生産を行う
        calc_list(a)        #　農家の生産状況を記録する
        Match(a, Demand)    #　均衡価格と需要量が決められる。
    if (t >= 11) and (t <= 20):
        tlist.append(t)     #　tをtlistに追加
        FarmerDecision(a)   #　農家が生産を行う

        calc_list(a)        #　農家の生産状況を記録する
        Match(a, Demand)    #　均衡価格と需要量が決められる。

#　シミュレーションパートの終わり

#　重要な変数のリストをそのまま表示
print("収入",profitlistA,"\n","\n","A生産量",qalistA,"\n","\n","自家労働力",LAlistA,"\n","\n","借りた労働力",Ln_relistA,"\n","\n",
      "貸した労働力",Ln_gilistA,"\n","\n", "余暇", llistA, "\n","\n", "自家農地", AAlistA, "\n","\n","借りた農地",An_relistA, "\n","\n",
      "貸した農地",An_gilistA, "\n","\n", "効用",utilitylistA, "\n","\n", "o消費",colistA,"\n","\n","C",ClistA, "\n","\n",
      "pA",pAlist,"\n","\n", "sellerの数",SellerNum,"\n","\n","buyerの数",BuyerNum)

#　sell&buy.py

4, 結果

結果は集計値のリストとして返されます。
めぼしいところとしては、総生産量、各種投入労働力、余暇、投入農地面積、効用、価格、net-sellerとbuyerの数、等です。

ちょっと全部考察すると面倒くさいので、価格と生産量と農家の種類だけ見てみます。

価格

↑matplotlibで簡単に出力した価格です。
実測値の初期値が0.6 peso/kg くらいなので、この条件で農家に生産計画を立てさせると、1.6ぐらいに価格が跳ね上がるということが分かりました。

ちなみに、ここにおける条件は、
①農家は4人いる
②一人は絶対にnet-sellerになるが、残り3人は最適な形態を選ぶ。
③保有農地面積にグラデーションがある（20ha, 2ha, 20ha, 200ha）
です。

あとは、価格が0になるときがありますよね。
これは、最適解の探索に失敗して、あり得ない解が選択され、使用可能な投入要素量を超えて生産計画が立てられ、生産量が爆増し、価格が落っこちた結果です。
ちなみに、こういうのにチクチク対応した結果、まあまあ落ち着いたのがこの結果です、、。

農家の最適化問題を解く部分を見てもらえたら、何かよくわかんない工程が挟まってますが、その理由は、変な値が出た時に、初期値を変えて探索し直させてるからです。

生産量

↑これは総生産量です。

ジグザグしてますね。この理由は、短期的に価格が増減を繰り返すので、それに反応して農家が生産量を調節したからです。
最後に爆増してるのは、価格のところで説明しましたが、最適化が失敗したからです。

今改めて見返してみると、変だな、、、。
需要の基準値が70,000kgなので、300,000kg作ったら、価格がかなり暴落するはずなのに、価格は微減でとどまってますね。何か悪いことが起きたのか、、、？

農家の種類

 sellerの数 [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2] 

 buyerの数 [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]

↑これは、農家がnet-seller（売る人）と、net-buyer（一部自給はするけど、売らない人）の数です。
農家は全部で4人いて、1人はnet-seller確定なので、もう一人だけ、最適解としてnet-sellerを選んだみたいですね。
で、生産者が不足した結果として価格が倍になったけど、それでもnet-buyerの人たちは農業が割に合わないと思ったみたいですね。

5, まとめ

ということで、今回はミクロな農家の生産計画と市場における価格決定を連関させた結果をお見せしました。
次の課題は、
①最適化問題の解き方を工夫する（詳しい人を見つけて、アドバイスをもらう）
②強化学習的なフレームも作ってみる
③ミクロな主体間の関連を増やす。例えば、雇われる人と雇う人をつなげるとか。
④マップを作って、その上に農地の区画を作って、どのエージェントがどの区画（農地）を所有しているか決めて、市場までの距離や雇い手・雇われる人の距離を決める。

です。
とりあえず、②に手を出してみようかな、、、。③と④は昔から温めてたことなので、こっちも捨てがたいですが、、、。
また、後で見返していて、たくさんコード的なミスや、私の不注意から生じた設定ミスを発見しました。次回はその辺も修正します。