#　はじめに

　近年、貧富の格差の拡大や富の分配、あるいは成長と分配といった言葉を良く耳にします。私たちの暮らしにおいても、消費税や所得税など税の負担は年々増していくように感じられます。国債残高が１千兆円を超え、支出と収入のバランスにも大きく課題を抱える日本経済ですが、MMT理論などの新しい貨幣理論も提唱されています。このような複雑で重要な経済を解き明かす学問が経済学ですが、プログラミングとシミュレーションによってより良く理解していきたいというのがこの記事の趣旨です。

　今回は、経済学の基本である需要と供給に関して、あるルールに基づいたシミュレーションにより、いわゆる「神の見えざる手」がいかに実現されるかを見てみたいと思います。そして、消費税が総生産や消費者余剰、生産者余剰に与える影響を考察してみたいと思います。

　また、今回のシミュレーションをインタラクティブに行えるWEBアプリを以下で公開しています。

Demand and supply balance: Streamlit

ミクロ経済のルール化

　ある市場で取引される財の量と価格がどのように決まるかを考えていきます。
まず、基本となる記号を定義していきましょう。

$P$: 価格
$X$: 供給量
$P=S(X)$: 供給曲線
$P=D(X)$: 需要曲線

　すなわち、供給曲線も需要曲線も量から価格を与える全単射の関数とします。これらからまずベースとなるルールAを定めます。

#####[ルールA]

供給量は供給曲線に従って決まる。すなわち、$X=S^{-1}(P)$。
価格は需要曲線に従って決まる。すなわち、$P=D(X)$。
供給量と価格は相互に影響しながら、収束するまで変化し続ける。

　1で供給量が供給曲線によって定められるのは、供給量を変更できる主体は生産者であり、供給曲線を知っているのは生産者のみであるためです。同様に、2で価格が需要曲線によって定められるのは、購買を決定するのは消費者であり、需要曲線を知っているのは消費者のみであるためです。

　ミクロ経済学の教科書の教えるところによれば、1～3の繰り返しにより、総余剰（消費者余剰＋生産者余剰）が最大となるという意味で最適な点、すなわち供給曲線と需要曲線が交わる点に収束するものとされています。いわゆる、アダム・スミスの言うところの「神の見えざる手」です。

　価格$P$と供給量$X$を時間依存の変数として、$P(n)$と$X(n)$を考えます。ただし、簡単のため$n$は離散時刻としましょう。このルールAを次に示すルールBに緩和します。

#####[ルールB]
1-a. 供給量は供給曲線に近づくよう変化する。
すなわち、

X(n+1) = \alpha X(n) + (1-\alpha) S^{-1}(P(n)), \\
1 \geq \alpha \geq 0

2-a. 価格は需要曲線に近づくよう変化する。
すなわち、

P(n+1) = \beta P(n) + (1-\beta) D(X(n)), \\
1 \geq \beta \geq 0

供給量と価格は相互に影響しながら、収束するまで変化し続ける。

　ここで、$\alpha,\beta$はモーメントの効果を表しています。以下のシミュレーションでは$\alpha=\beta=0.9$を用います。

　さて、上記の計算を行うには供給曲線と需要曲線の具体形を与える必要があります。ここでは、以下の式を仮定しましょう。

\begin{align}
S(X) &= S_0 exp(s X), s>0, \\
D(X) &= D_0 exp(d X), d<0, \\
D_0 &> S_0
\end{align}

このとき、$S^{-1}$は次式となります。

S^{-1}(P) = ln(\frac{P}{S_0})/s

#税効果を入れてみる

　さらに、消費税の効果を入れてみましょう。このとき、ルールBに税効果を入れたルールCを考えます。

#####[ルールC]
1-b. 供給量は市場価格$P(n)$から消費税率$r_C$を割り引いた価格に基づいた供給曲線に近づくよう変化する。すなわち、

X(n+1) = \alpha X(n) + (1-\alpha) S^{-1}(\frac{P(n)}{1+r_C}), \\
1 \geq \alpha \geq 0

2-b. 2-aと同じ
3. 同上

総余剰を計算する

　ある時点nで市場価格が$P(n)$、供給量が$X(n)$であるとき、消費者余剰$CS(n)$、生産者余剰$PS(n)$および総余剰$TS(n)$を次式で定義します。消費者余剰とは、消費者にとって市場と需要曲線のギャップから受ける利益、生産者余剰とは、生産者にとって市場と供給曲線のギャップから受ける利益、総余剰はそれらの和と言えます。

\begin{align}
CS(n) &= \int_0^{X(n)} (D(x)-P(n)) dx, \\
PS(n) &= \int_0^{X(n)} (P(n)-S(x)) dx, \\
TS(n) &= CS(n)+PS(n)
\end{align}

具体的に計算すると、

\begin{align}
CS(n) &= \int_0^{X(n)} (D(x)-P(n)) dx, \\
      &= \int_0^{X(n)} (D_0 exp(d \cdot x) - P(n)) dx, \\
      &= \frac{D_0}{d} (exp(d X(n)) - 1) - P(n) X(n), \\
PS(n) &= \int_0^{X(n)} (P(n)-S(x)) dx, \\
      &= \int_0^{X(n)} (P(n) - S_0 exp(s \cdot x) ) dx, \\
      &= P(n) X(n) - \frac{S_0}{s}(exp(s X(n)) - 1), \\
TS(n) &= CS(n)+PS(n), \\
      &= \frac{D_0}{d} (exp(d X(n)) - 1) - \frac{S_0}{s} (exp(s X(n)) - 1)
\end{align}

となります。

Pythonコード

　以上の計算をPythonで実装してみましょう。まずは、計算パートです。

############################################
# class
class MicroModel01:
    def __init__(self):
        # supply curve
        self.S_0 = 20.0
        self.s = 0.016
        self.alpha = 0.9
        # demand curve
        self.D_0 = 100.0
        self.d = -0.025
        self.beta = 0.9
        # tax
        self.r_C = 0.1
    
    def surplus(self, P, X):
        CS = self.D_0 / self.d * (np.exp(self.d * X) - 1) - P * X
        PS = P * X - self.S_0 / self.s * (np.exp(self.s * X) - 1)
        TS = CS+PS
        return CS, PS, TS
    
    def demand(self, X):
        return self.D_0 * np.exp( self.d * X)
    
    def supply(self, X):
        return self.S_0 * np.exp( self.s * X)
    
    def invsupply(self, P):
        return np.log(P/self.S_0)/self.s

############################################
# sim functions
def calcsim(m, N, X0, P0):
    #
    hn = [0]
    hx = [X0]
    hp = [P0]
    x = X0
    p = P0
    #
    cs, ps, ts = m.surplus(p,x)
    hcs = [cs]
    hps = [ps]
    hts = [ts]
    #
    for n in range(1,N):
        # 1. supply
        xn = m.alpha * x + (1- m.alpha) * m.invsupply( p / (1 + m.r_C))
        # 2. demand
        pn = m.beta * p + (1 - m.beta) * m.demand( x )
        #
        x = xn
        p = pn
        cs, ps, ts = m.surplus(p,x)
        #
        hn.append(n)
        hx.append(x)
        hp.append(p)
        hcs.append(cs)
        hps.append(ps)
        hts.append(ts)
    return hn, hx, hp, hcs, hps, hts

def taxsim():
    N= 100
    X0 = 20
    P0 = 80
    m = MicroModel01()
    #
    hr = []
    hfx = []
    hfp = []
    hfcs = []
    hfps = []
    hfts = []
    hgdp = []
    hinc = []
    htax = []
    #
    for r in range(0,30):
        m.r_C = 0.01*r
        hn, hx, hp, hcs, hps, hts = calcsim(m, N, X0, P0)
        hr.append(r)
        hfx.append(hx[-1])
        hfp.append(hp[-1])
        hfcs.append(hcs[-1])
        hfps.append(hps[-1])
        hfts.append(hts[-1])
        hgdp.append(hx[-1]*hp[-1])
        hinc.append(hx[-1]*hp[-1]*(1-m.r_C))
        htax.append(hx[-1]*hp[-1]*m.r_C)
    #
    return hr, hfx, hfp, hfcs, hfps, hfts, hgdp, hinc, htax

次に、グラフ化パートです。

############################################
# draw functions
def drawcurve(f, xmin, xmax, label=None):
    N = 100
    x = [(xmax - xmin)*i/(N-1) + xmin for i in range(N)]
    y = [f(v) for v in x]
    plt.plot(x,y,label=label)

def drawtrend(m, hn, hx, hp, hcs, hps, hts):
    fig1 = plt.figure(figsize=(6,7))
    plt.subplot(211)
    plt.plot(hn ,hx, label='X')
    plt.plot(hn ,hp, label='P')
    plt.xlabel("time")
    plt.legend()
    plt.subplot(212)
    plt.plot(hn, hcs, label='CS')
    plt.plot(hn, hps, label='PS')
    plt.plot(hn, hts, label='TS')
    plt.xlabel("time")
    plt.legend()

def drawhis(m, hn, hx, hp):
    fig = plt.figure()
    drawcurve(m.demand, 0, 100)
    drawcurve(m.supply, 0, 100)
    plt.plot(hx ,hp, label='sim')
    plt.xlabel("amount")
    plt.ylabel("price")
    plt.legend()

def drawtaxsim(hr, hfx, hfp, hfcs, hfps, hfts, hgdp, hinc, htax):
    fig1 = plt.figure(figsize=(6,12))
    plt.subplot(411)
    plt.plot(hr, hfx, label='X')
    plt.plot(hr, hfp, label='P')
    plt.xlabel("consumption tax rate [%]")
    plt.legend()
    plt.subplot(412)
    plt.plot(hr, hgdp, label='Gross Production')
    plt.legend()
    plt.subplot(413)
    plt.plot(hr, hgdp, label='GP')
    plt.plot(hr, hinc, label='Income')
    plt.plot(hr, htax, label='Tax')
    plt.legend()
    plt.subplot(414)
    plt.plot(hr, hfcs, label='CS')
    plt.plot(hr, hfps, label='PS')
    plt.plot(hr, hfts, label='TS')
    plt.legend()

シミュレーション

　以上の準備のもと、いくつかのケースでシミュレーションを実行してみましょう。

ケース１ (消費税がない場合)

　まず、消費税がない場合$(r_C= 0 )$の平衡に至るまでの過程をシミュレーションします。初期点は適当に$(X_0, P_0)=(20, 80)$へと置いておきます。

m0 = MicroModel01()
m0.r_C = 0
X0, P0 = 20, 80
hn, hx, hp, hcs, hps, hts = calcsim(m0, 100, X0, P0)
drawtrend(m0, hn, hx, hp, hcs, hps, hts)
drawhis(m0, hn, hx, hp)

　上記のグラフから次の傾向が見て取れます。

収束の過程では時計回りの曲線を描く。
需要曲線と供給曲線の交点に収束し、総余剰が最大値となる。

ケース２ (消費税が10%の場合)

　次に、消費税が10%の場合$(r_C= 0.1 )$の平衡に至るまでの過程をシミュレーションします。初期点は$(X_0, P_0)=(20, 80)$と置いておきます。

m1 = MicroModel01()
X0, P0 = 20, 80
hn, hx, hp, hcs, hps, hts = calcsim(m1, 100, X0, P0)
drawtrend(m1, hn, hx, hp, hcs, hps, hts)
drawhis(m1, hn, hx, hp)

　ケース１と比較すると、次の傾向が見て取れます。

時計回りの形状は同じ。
収束点が需要曲線と供給曲線の交点から需要曲線上の左上にずれる。
総余剰が減少する。

ケース３ (初期点が異なる場合)

　消費税は10%のまま、初期点を$(X_0, P_0)=(80, 80)$に変えてみます。

X0, P0 = 80, 80
hn, hx, hp, hcs, hps, hts = calcsim(m1, 100, X0, P0)
drawtrend(m1, hn, hx, hp, hcs, hps, hts)
drawhis(m1, hn, hx, hp)

　ケース２と比較すると、次の傾向が見て取れます。

収束点は同じであるものの、途中の過程でCSが大きく下がり負の値となる。また、PSも一旦0付近まで下がる。

　つまり、初期点によって収束点は変わらないものの、途中の過程が変わり、CSやPSが一時的に大きく変化することがあることが分かります。

ケース４消費税の影響)

　消費税率を0%から30%まで変化させて、収束点がどのように変化するのか見てみましょう。

hr, hfx, hfp, hfcs, hfps, hfts, hgdp, hinc, htax = taxsim()
drawtaxsim(hr, hfx, hfp, hfcs, hfps, hfts, hgdp, hinc, htax)

　このグラフから、次の傾向が見て取れます。

消費税率が上昇すると、収束時の価格が上昇し、供給量が減少する。
消費税率が上昇すると、総生産(価格×供給量)は低下するが、消費税収入は税率にほぼ比例して上昇する。
消費税率が上昇すると、消費者余剰（CS)は低下するが、生産者余剰（PS)は逆に上昇する。

まとめ

　以上のシミュレーションから、ミクロ経済学シミュレーションに関して以下の傾向が分かりました。

供給量は供給曲線に、価格は需要曲線にそれぞれ近づくという過程の繰り返しにより、総余剰が最大となる、供給曲線と需要曲線の交点に任意の初期点から収束していく。
神の見えざる手（収束の軌跡）は時計回りである。
初期点(X0, P0)に関して、P0が同一の場合、供給少の場合の方が、供給大の場合よりも、収束過程の消費者余剰（CS)と生産者余剰（PS)が高く保たれる場合がある。
消費税率が上昇すると、総生産は低下し、消費者余剰（CS)は低下する。一方、生産者余剰は上昇する場合がある。

参考

　下記の文献等を参考にしました。

ミクロ経済学の力
 アダム･スミスが経済学を語る上で外せない訳