RBCモデルによる持ち家率の簡素な分析

RBCモデル

• Real Business Cycle = 実物的景気循環
• 技術進歩、生産性、嗜好などの実物的ショックが景気循環の要因
• 家計・企業は合理的期待のもとで最適に行動（効用最大化・利潤最大化）
  • 今回は企業を明示的に数理化してないので厳密にはRBCモデルではない
• 価格・賃金は完全に柔軟
• 更に経済事象を組み込んだDSGEモデル(Dynamic Stochastic General Equilibrium Model = 動学的確率的一般均衡モデル)は各国中央銀行や学界などが開発し利用

家計(効用最大化)

\begin{align}
\max_{\{c_t, l_t, h_t, \theta_t, k_{t+1}, h_{t+1}\}_{t=0}^{\infty}} 
& \quad 
\mathbb{E}_t \sum_{t=0}^{\infty} \beta^t \, \xi_t 
\left[ 
\frac{c_t^{1-\sigma}}{1-\sigma} 
- \frac{\mu \, l_t^{1+\lambda}}{1+\lambda} 
+ \psi \frac{(\theta_t \, h_t)^{1-\eta}}{1-\eta} 
\right] 
\\[2mm]
\text{s.t.} \quad
& c_t + q_t \, (h_{t+1} - (1-\delta_h) h_t) + k_{t+1} - (1-\delta_k) k_t 
= w_t \, l_t + R_t \, (1-\theta_t) \, h_t + r_t \, k_t, \quad
\\[1mm]
& v_t k_t^\alpha l_t^{1-\alpha} = w_t l_t + r_t k_t \quad
\\[1mm]
& k_{t+1} = (1-\delta_k) k_t + i_k \quad  \\
& h_{t+1} = (1-\delta_h) h_t + i_h \quad 
\end{align}

ジャンプ変数

• $c_t$ ：消費
• $l_t$ ：労働供給
• $q_t$ ：住宅価格
• $\theta_t$ ：持ち家率（自己使用住宅の比率）
• $R_t$ ：賃貸住宅サービス価格

状態変数

• $k_t$ ：物的資本ストック
• $h_t$ ：住宅ストック
• $v_t$ ：技術ショック
• $\xi_t$ ：嗜好ショック（効用の重み）

外生変数

• $r_t$ :物的資本のレンタル価格
• $w_t$ :賃金

ショック

• $\varepsilon^v_t$ ：技術ショックのイノベーション
• $\varepsilon^\xi_t$ ：嗜好ショックのイノベーション

パラメータ

• $\beta \in (0,1)$ ：割引因子
• $\sigma > 0$ ：相対的リスク回避度
• $\mu > 0$ ：労働の不効用の重み
• $\lambda > 0$ ：逆フリッシュ弾力性
• $\psi > 0$ ：住宅サービス効用ウェイト
• $\eta > 0$ ：住宅サービス効用の曲率
• $\alpha \in (0,1)$ ：資本分配率
• $\delta_k \in (0,1)$ ：物的資本の減耗率
• $\delta_h \in (0,1)$ ：住宅の減耗率
• $\phi \in (0,1)$ ：技術ショックの自己相関
• $\rho_\xi \in (0,1)$ ：嗜好ショックの自己相関
• $\theta^{ss} \in (0,1)$ ：定常状態の持ち家率
• $\chi_q$ ：住宅価格感応度
• $\chi_R$ ：賃貸価格感応度
• $\chi_\xi$ ：嗜好ショック感応度

最適化問題の解法(ラグランジュ関数を定義して紙と鉛筆で解く)

\begin{align}
\mathcal{L} = & \, \mathbb{E}_t \sum_{t=0}^{\infty} \beta^t \, \xi_t \left[ \frac{c_t^{1-\sigma}}{1-\sigma} - \frac{\mu l_t^{1+\lambda}}{1+\lambda} + \psi \frac{(\theta_t h_t)^{1-\eta}}{1-\eta} \right] \\
& + \sum_{t=0}^{\infty} \lambda_t \left[ w_t l_t + R_t (1-\theta_t) h_t + r_t k_t - c_t - q_t \left( h_{t+1} - (1-\delta_h) h_t \right) - k_{t+1} + (1-\delta_k) k_t \right] \\
& + \sum_{t=0}^{\infty} \mu_t \left[ k_{t+1} - (1-\delta_k) k_t - i_k \right] \\
& + \sum_{t=0}^{\infty} \nu_t \left[ h_{t+1} - (1-\delta_h) h_t - i_h \right] \\
& + \sum_{t=0}^{\infty} \iota_t \left[ v_t k_t^\alpha l_t^{1-\alpha} -  w_t l_t - r_t k_t \right]
\end{align}

非線形経済システム = First Order Condition(1階条件) + アドホックルール

\xi_t c_t^{-\sigma}
=
\beta E_t \Big[
\xi_{t+1} c_{t+1}^{-\sigma}
\big(
1-\delta_k
+ \alpha v_{t+1} k_{t+1}^{\alpha-1} l_{t+1}^{1-\alpha}
\big)
\Big]

\mu l_t^{\lambda}
=
\xi_t c_t^{-\sigma}
(1-\alpha) v_t k_t^{\alpha} l_t^{-\alpha}

q_t
=
\beta E_t
\left[
\frac{\xi_{t+1} c_{t+1}^{-\sigma}}{\xi_t c_t^{-\sigma}}
\left(
(1-\delta_h) q_{t+1}
+
\psi \xi_{t+1} (\theta_{t+1} h_{t+1})^{-\eta}
\right)
\right]

\psi \xi_t (\theta_t h_t)^{-\eta}
=
q_t c_t^{-\sigma}

v_t k_t^{\alpha} l_t^{1-\alpha}
=
c_t
+ q_t \bigl(h_{t+1} - (1-\delta_h) h_t\bigr)
+ k_{t+1}
- (1-\delta_k) k_t

\theta_t - \theta^{ss}
=
\chi_q (q_t - q^{ss})
- \chi_R (R_t - R^{ss})
- \chi_\xi (\xi_t - \xi^{ss})

v_{t+1} = \phi v_t + \varepsilon^v_t

\xi_{t+1} = \rho_{\xi} \xi_t + \varepsilon^{\xi}_t

コード

python

# -*- coding: utf-8 -*-
import numpy as np
from scipy.linalg import ordqz
import matplotlib.pyplot as plt

# -------------------------------
# Numerical Jacobian (central differences)
# -------------------------------
def numerical_jacobian(func, x0, eps=1e-6):
    """
    Compute a central-difference Jacobian J of func at x0.
    func: R^n -> R^m, returns a vector residual
    J: (m x n) matrix of partial derivatives
    """
    x0 = np.asarray(x0, float)
    f0 = np.asarray(func(x0))
    m, n = f0.size, x0.size
    J = np.zeros((m, n))
    h = eps * (1 + np.abs(x0))
    for j in range(n):
        x1 = x0.copy();  x1[j] += h[j]
        x2 = x0.copy();  x2[j] -= h[j]
        f1 = np.asarray(func(x1))
        f2 = np.asarray(func(x2))
        J[:, j] = (f1 - f2) / (2 * h[j])
    return J

# -------------------------------
# Residual builder: bind names -> vector-valued function
# -------------------------------
def make_residual_function(varnames, shock_names, residual_eqs):
    """
    Build a residual function F(x, params) from named variables and equations.
    x packs [future vars, current vars, shocks]
    """
    def residual(x, params):
        X = {name: x[i] for i, name in enumerate(varnames + shock_names)}
        return np.array([eq(X, params) for eq in residual_eqs])
    return residual

# -------------------------------
# QZ solver with BK check (ordqz(C, B))
# -------------------------------
def solve_qz_system(B, C, njump, npred, D=None, verbose=False):
    """
    Solve B x_{t+1} + C x_t + D eps_t = 0 via generalized Schur (QZ).
    Ordering uses ordqz(C, B) so eigenvalues are lambda = alpha/beta of (C - lambda B).
    BK condition: number of unstable roots must equal njump.
    Returns policy P (jump = P state), transition A (state), shock map G.
    """

    def sort_unstable(alpha, beta):
        # Mark eigenvalues outside unit circle as "unstable"
        tol = 1e-12
        alpha = np.asarray(alpha)
        beta  = np.asarray(beta)
        mask_zero_beta = np.abs(beta) < tol
        lam_abs = np.zeros_like(alpha, float)
        lam_abs[mask_zero_beta] = np.inf
        lam_abs[~mask_zero_beta] = np.abs(alpha[~mask_zero_beta] / beta[~mask_zero_beta])
        return lam_abs > 1.0

    # Run QZ with reordering by custom callable
    res = ordqz(C, B, sort=sort_unstable, output="real")

    # Handle SciPy return formats (object or tuple)
    if hasattr(res, "AA"):
        AA = res.AA; BB = res.BB; alpha = res.alpha; beta = res.beta; Q = res.Q; Z = res.Z
    else:
        seq = list(res)
        if len(seq) == 6:
            AA, BB, alpha, beta, Q, Z = seq
        elif len(seq) == 4:
            AA, BB, Q, Z = seq
            alpha = np.diag(AA); beta = np.diag(BB)
        else:
            raise RuntimeError(f"Unexpected ordqz output length {len(seq)}")

    lam = np.array(alpha) / np.array(beta)
    mask_unstable = np.abs(lam) > 1.0
    nunstable = int(mask_unstable.sum())
    if verbose:
        print("QZ generalized eigenvalues:", lam)
        print("Unstable roots:", nunstable, "Jump vars:", njump)
    if nunstable != njump:
        raise RuntimeError(f"Blanchard-Kahn failed: unstable roots={nunstable} != jump vars={njump}")

    # Stable invariant subspace (right Schur vectors)
    Z_stable = Z[:, ~mask_unstable]

    # Variable ordering: [jump | state] after grouping
    jump_idx  = list(range(njump))             # first njump rows
    state_idx = list(range(njump, njump+npred))# next npred rows

    # Policy: jump = P * state
    ZJ = Z_stable[jump_idx, :]
    ZS = Z_stable[state_idx, :]
    P = (np.linalg.solve(ZS.T, ZJ.T)).T

    # State transition A on the stable block
    idx_stable = np.where(~mask_unstable)[0]
    AA22 = AA[np.ix_(idx_stable, idx_stable)]
    BB22 = BB[np.ix_(idx_stable, idx_stable)]
    A_schur = np.linalg.solve(BB22, AA22)
    A = ZS @ A_schur @ np.linalg.inv(ZS)

    # Shock map G
    if D is None:
        G = np.zeros((npred, 1))
    else:
        D_q = Q.T @ D          # left transform to QZ coordinates
        D2 = D_q[idx_stable, :]# keep stable rows
        G_schur = np.linalg.solve(BB22, D2)
        G = ZS @ G_schur

    return P.real, A.real, G.real, lam.real

# -------------------------------
# Linear simulation: state and jump paths
# -------------------------------
def simulate_system(A, G, P, T, shock_process, x_state0=None):
    """
    Simulate linear system:
      state_{t+1} = A state_t + G eps_t
      jump_t      = P state_t
    """
    nstate = A.shape[0]
    njump  = P.shape[0]
    x_state = np.zeros(nstate) if x_state0 is None else np.asarray(x_state0, float)

    state_path = np.zeros((T+1, nstate))
    jump_path  = np.zeros((T+1, njump))
    state_path[0] = x_state
    jump_path[0]  = P @ x_state

    for t in range(T):
        eps = shock_process(t) if shock_process is not None else np.zeros(G.shape[1])
        x_next = A @ x_state + G @ eps
        state_path[t+1] = x_next
        jump_path[t+1]  = P @ x_next
        x_state = x_next

    return state_path, jump_path

# -------------------------------
# Build model: steady state, residuals, Jacobian, blocks
# -------------------------------

def build_model(params):
    def ss_vals(p):
        v = 1.0; xi = 1.0; theta = p["theta_ss"]
        kls = (((1/p["beta"]) + p["delta_k"] - 1) / p["alpha"])**(1/(p["alpha"]-1))
        w = (1-p["alpha"]) * kls**p["alpha"]
        cl = kls**p["alpha"] - p["delta_k"]*kls
        l = ((w/p["mu"]) * (cl**(-p["sigma"])) * xi)**(1/(p["lamda"]+p["sigma"]))
        k = l*kls
        c = l*cl
        h = ((p["psi"]*(1-p["beta"]))/(p["beta"]*(1-p["delta_h"])))**(1/p["eta"]) * c**(p["sigma"]/p["eta"])
        q = (p["beta"] * p["psi"] * xi * (theta*h)**(-p["eta"])) / (1 - p["beta"] * (1 - p["delta_h"]))
        R = (p["psi"] * xi * ((1.0-theta)*h)**(-p["eta"])) / (c**(-p["sigma"]))
        return dict(l=l, c=c, k=k, h=h, q=q, R=R, v=v, xi=xi, theta=theta)

    ss = ss_vals(params)

    # --- REORDER to [jump | state] in BOTH future & current ---
    future_vars  = ["l_f","c_f","q_f","theta_f","k_f","h_f","v_f","xi_f"]
    current_vars = ["l","c","q","theta","k","h","v","xi"]
    shock_vars   = ["eps_v","eps_xi"]
    varnames = future_vars + current_vars

    # Equations (unchanged except theta_f is used in price Euler)
    def eq1(X,p):
        mpk_f = p["alpha"]*X["k_f"]**(p["alpha"]-1) * X["l_f"]**(1-p["alpha"]) * X["v_f"]
        return X["xi"]*X["c"]**(-p["sigma"]) - p["beta"]*X["xi_f"]*X["c_f"]**(-p["sigma"]) * (1-p["delta_k"] + mpk_f)

    def eq2(X,p):
        w_t = (1-p["alpha"]) * X["k"]**p["alpha"] * X["l"]**(-p["alpha"]) * X["v"]
        return -p["mu"]*X["l"]**p["lamda"] + X["xi"]*X["c"]**(-p["sigma"]) * w_t

    def eq3(X,p):
        y = X["k"]**p["alpha"] * X["l"]**(1-p["alpha"]) * X["v"]
        return y - X["c"] - X["q"]*(X["h_f"] - (1-p["delta_h"]) * X["h"]) - X["k_f"] + (1-p["delta_k"]) * X["k"]

    def eq4(X,p):
        m = (X["xi_f"]*X["c_f"]**(-p["sigma"]))/(X["xi"]*X["c"]**(-p["sigma"]))
        return X["q"] - p["beta"]*m*((1-p["delta_h"])*X["q_f"] + p["psi"]*X["xi_f"]*(X["theta_f"]*X["h_f"])**(-p["eta"]))

    def eq5(X,p):
        return X["v_f"] - p["phi"] * X["v"] - X["eps_v"]

    def eq6(X,p):
        return X["xi_f"] - p["rho_xi"] * X["xi"] - X["eps_xi"]

    def eq7(X,p):
        return p["psi"]*X["xi"]*(X["theta"]*X["h"])**(-p["eta"]) - X["q"]*X["c"]**(-p["sigma"])

    def eq8(X,p):
        R_cur = (p["psi"]*X["xi"]*((1.0-X["theta"])*X["h"])**(-p["eta"]))/(X["c"]**(-p["sigma"]))
        return (X["theta"] - p["theta_ss"]) - p["chi_q"]*(X["q"] - ss["q"]) + p["chi_R"]*(R_cur - ss["R"]) + p["chi_xi"]*(X["xi"] - ss["xi"])

    residual_eqs = [eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8]
    residual = make_residual_function(varnames, shock_vars, residual_eqs)

    # Steady-state vector consistent with new ordering
    xss = np.array([
        ss["l"], ss["c"], ss["q"], ss["theta"], ss["k"], ss["h"], ss["v"], ss["xi"],
        ss["l"], ss["c"], ss["q"], ss["theta"], ss["k"], ss["h"], ss["v"], ss["xi"],
        0.0, 0.0
    ], float)

    # Jacobian
    J = numerical_jacobian(lambda x: residual(x, params), xss)
    J_fc = J[:, :16]

    # Column scaling (future/current each length-8, same ordering)
    scale_f = np.array([ss["l"], ss["c"], ss["q"], ss["theta"], ss["k"], ss["h"], ss["v"], ss["xi"]])
    scale_c = np.array([ss["l"], ss["c"], ss["q"], ss["theta"], ss["k"], ss["h"], ss["v"], ss["xi"]])
    TW_f = np.tile(scale_f, (len(residual_eqs), 1))
    TW_c = np.tile(scale_c, (len(residual_eqs), 1))

    B = -J_fc[:, :8]   * TW_f
    C =  J_fc[:,  8:16]* TW_c

    D_raw = J[:, -2:]
    D = D_raw * np.array([ss["v"], ss["xi"]])

    return B, C, D, ss

# -------------------------------
# Top-level solve wrapper
# -------------------------------
def try_solve(params, verbose=False):
    """
    Build model and solve QZ with BK check. Expect njump=4 (l,c,q,theta), npred=4 (k,h,v,xi).
    """
    B, C, D, ss = build_model(params)
    if verbose:
        print("Shapes -> B:", B.shape, "C:", C.shape, "D:", D.shape)
    P, A, G, lam = solve_qz_system(B, C, njump=4, npred=4, D=D, verbose=verbose)
    return (P, A, G, lam), ss

# -------------------------------
# Helpers: rent series and plotting
# -------------------------------
def rent_series(xi_path, c_path, theta_path, h_path, psi, eta, sigma):
    """
    Rental service price implied by FOC:
      R_t = [psi * xi_t * ((1 - theta_t) h_t)^(-eta)] / [c_t^(-sigma)]
    """
    return (psi*xi_path*((1.0-theta_path)*h_path)**(-eta)) / (np.maximum(1e-12, c_path**(-sigma)))


def save_irfs(prefix, state_path, jump_path, ss, params):
    """
    Plot IRFs using level-safe reconstruction:
    - Rebuild levels from hat variables: X_level = X_ss * (1 + X_hat)
    - Clamp theta in (0,1) and h positive to avoid NaNs in R
    - Plot rent as percentage deviation: R_hat = R_level / R_ss - 1
    """

    # Unpack hat (deviation) series from paths
    k_hat  = state_path[:,0]
    h_hat  = state_path[:,1]
    v_hat  = state_path[:,2]
    xi_hat = state_path[:,3]

    l_hat     = jump_path[:,0]
    c_hat     = jump_path[:,1]
    q_hat     = jump_path[:,2]
    theta_hat = jump_path[:,3]

    # Rebuild levels from hat variables: X_level = X_ss * (1 + X_hat)
    k_lev   = ss["k"]   * (1.0 + k_hat)
    h_lev   = ss["h"]   * (1.0 + h_hat)
    v_lev   = ss["v"]   * (1.0 + v_hat)
    xi_lev  = ss["xi"]  * (1.0 + xi_hat)
    l_lev   = ss["l"]   * (1.0 + l_hat)
    c_lev   = ss["c"]   * (1.0 + c_hat)
    q_lev   = ss["q"]   * (1.0 + q_hat)
    # Theta must remain in (0,1): clip to avoid ((1-theta)*h)^(-eta) NaNs
    theta_lev = np.clip(ss["theta"] * (1.0 + theta_hat), 1e-8, 1.0 - 1e-8)

    # Level-safe rent from FOC
    psi   = params["psi"]
    eta   = params["eta"]
    sigma = params["sigma"]
    R_lev = (psi * xi_lev * ((1.0 - theta_lev) * h_lev)**(-eta)) / (np.maximum(1e-12, c_lev**(-sigma)))
    # Rent IRF as percentage deviation
    R_hat = (R_lev / ss["R"]) - 1.0

    # Optional: console diagnostics to ensure we are not plotting NaNs
    print(f"[{prefix}] R_hat min={np.nanmin(R_hat):.6f}, max={np.nanmax(R_hat):.6f}")

    # For consistency, plot other series also as percentage deviations
    k_hat_pct  = (k_lev  / ss["k"])  - 1.0
    h_hat_pct  = (h_lev  / ss["h"])  - 1.0
    q_hat_pct  = (q_lev  / ss["q"])  - 1.0
    c_hat_pct  = (c_lev  / ss["c"])  - 1.0
    l_hat_pct  = (l_lev  / ss["l"])  - 1.0
    theta_hat_pct = (theta_lev / ss["theta"]) - 1.0  # “rate” shown as deviation from ss

    t = np.arange(len(k_hat))

    # Homeownership (as deviation)
    plt.figure(figsize=(6,4)); plt.plot(t, theta_hat_pct)
    plt.title(prefix+"_homeownership"); plt.grid(True); plt.tight_layout()
    plt.savefig(prefix+"_homeownership.png", dpi=150); plt.close()

    # Housing stock (deviation)
    plt.figure(figsize=(6,4)); plt.plot(t, h_hat_pct)
    plt.title(prefix+"_housing_stock"); plt.grid(True); plt.tight_layout()
    plt.savefig(prefix+"_housing_stock.png", dpi=150); plt.close()

    # House price q (deviation)
    plt.figure(figsize=(6,4)); plt.plot(t, q_hat_pct)
    plt.title(prefix+"_house_price_q"); plt.grid(True); plt.tight_layout()
    plt.savefig(prefix+"_house_price_q.png", dpi=150); plt.close()

    # Rent R (deviation)
    plt.figure(figsize=(6,4)); plt.plot(t, R_hat)
    plt.title(prefix+"_rent_R"); plt.grid(True); plt.tight_layout()
    plt.savefig(prefix+"_rent_R.png", dpi=150); plt.close()

    # Consumption (deviation)
    plt.figure(figsize=(6,4)); plt.plot(t, c_hat_pct)
    plt.title(prefix+"_consumption"); plt.grid(True); plt.tight_layout()
    plt.savefig(prefix+"_consumption.png", dpi=150); plt.close()

    # Labor (deviation)
    plt.figure(figsize=(6,4)); plt.plot(t, l_hat_pct)
    plt.title(prefix+"_labor"); plt.grid(True); plt.tight_layout()
    plt.savefig(prefix+"_labor.png", dpi=150); plt.close()

    # Capital (deviation)
    plt.figure(figsize=(6,4)); plt.plot(t, k_hat_pct)
    plt.title(prefix+"_capital"); plt.grid(True); plt.tight_layout()
    plt.savefig(prefix+"_capital.png", dpi=150); plt.close()

# -------------------------------
# IRF simulation (unit shock at t=0)
# -------------------------------
def simulate_irf(A, G, P, T, amp, shock):
    """
    Produce an IRF for either preference ('pref') or technology ('tech') shock.
    eps = [eps_v, eps_xi], with a unit shock in selected dimension at t=0.
    """
    if shock=="pref":
        def shock_process(t): return np.array([0.0, amp if t==0 else 0.0])
    else:
        def shock_process(t): return np.array([amp if t==0 else 0.0, 0.0])
    state_path, jump_path = simulate_system(A, G, P, T, shock_process)
    return state_path, jump_path

# -------------------------------
# Main
# -------------------------------
if __name__ == "__main__":
    print("Building tenure-augmented housing model: theta is endogenous (jump) and reacts to shocks.")
    params = dict(
        sigma=2.0, alpha=0.33, beta=0.99,
        delta_k=0.025, delta_h=0.02,
        mu=1.0, lamda=1.0, psi=0.06, eta=1.7,
        phi=0.85, rho_xi=0.85,
        theta_ss=0.6,
        chi_q=-0.10, chi_R=0.10, chi_xi=0.05
    )

    print("Constructing Jacobian and scaling: future=8, current=8; solving QZ (expect jumps=4, states=4) ...")
    try:
        (P, A, G, lam), ss = try_solve(params, verbose=True)
        print("BK OK. Shapes -> P", P.shape, "A", A.shape, "G", G.shape)
        print("Steady state:", ss)

        T = 40; amp = 0.01
        print("Simulating preference shock IRF and saving figures ...")
        sp, jp = simulate_irf(A, G, P, T, amp, "pref")
        save_irfs("irf_pref", sp, jp, ss, params)

        print("Simulating technology shock IRF and saving figures ...")
        st, jt = simulate_irf(A, G, P, T, amp, "tech")
        save_irfs("irf_tech", st, jt, ss, params)

        print("Done. Saved: irf_pref_*.png, irf_tech_*.png")

    except RuntimeError as e:
        print("BK failed. Tune parameters: reduce |chi_q| and chi_R, raise eta, lower phi/rho_xi slightly.")
        print("Error detail:", e)

処理手順

1 パラメータ設定と定常状態の計算
2 残差関数とヤコビアンの構築
3 QZ分解とBK条件に基づいた動学的特性の解析
4 ショックに対するインパルス応答関数（IRF）のシミュレーション
5 結果の保存と可視化

結果(インパルスレスポンス)

• 縦軸×100 = %乖離(例: 縦軸の値 0.01 → 定常状態から1%乖離と解釈)
  

結論

• 嗜好ショック→持ち家率: 最初は増えそれ以降減少
  • 住宅サービスが欲しい→住宅所有比率を上げる(持ち家率の上昇)→住宅投資が増える→住宅供給が増える→家賃が減る + 嗜好ショックが減る + 住宅価格が減る→ 持ち家率の低下

今後

• 賃金・資本レンタル率の内生化(利潤最大化する企業の追加)
• 家計に異質性を導入
• 国内総生産(=総供給)の表示
• 社会計画者の問題として設定
• 厚生分析
• 価格硬直性、独占的競争市場、金融市場の不完全性を導入
• 政府・中央銀行の導入
  • 財政・金融政策の効果検証
• Occasionally Binding Constraints(時々かかる制約)
  • 借入制約、ゼロ金利制約(ZLB)など
• Dynareの活用
  • 公式
    https://www.dynare.org/
  • マニュアル
    https://www.dynare.org/manual.pdf
  • Ubuntu on IBM CloudにOctaveとdynareをインストール + 経済モデルシミュレーション
    https://qiita.com/Yuta_Nakamura_/items/9ffb8c554483b89960ed　

参考

• 現代マクロ経済学講義―動学的一般均衡モデル入門 加藤 涼
  www.amazon.co.jp/dp/4492313702
• DSGEモデルによるマクロ実証分析の方法 廣瀬康生
  www.amazon.co.jp/dp/4943852394
• 動学的一般均衡モデルによる財政政策の分析 江口允崇
  www.amazon.co.jp/dp/4943852343