お勉強メモ

経済学・計量経済学・統計学などのお勉強メモです。

経済学メモ:ARモデルのインパルス応答関数

経済学 計量経済学 時系列分析 Python

◆本稿の内容
AR(1)モデルとAR(2)モデルのインパルス応答関数(以下ではIRFと表記します)を導出します。実データを用いて、AR(1)モデルとAR(2)モデルとそれぞれのIRFの推定も行います。

◆目次

Ⅰ:AR(1)モデルのIRF

以下のAR(1)モデルを考える。ただし、|a_{1}|<1、\varepsilon_{t}はホワイトノイズとする。


y_{t}=a_{0}+a_{1}y_{t-1}+\varepsilon_{t}

このAR(1)モデルを解くと、


y_{t}=\dfrac{a_{0}}{1-a_{1}}+\sum\limits_{i=0}^{\infty} a_{1}^{i}\varepsilon_{t-i}・・・①

となる。*1以下では、\mu=\dfrac{a_{0}}{1-a_{1}}と書く。①式を\varepsilon_{t}で微分すると、


\begin{eqnarray} \dfrac{\partial y_{t}}{\partial \varepsilon_{t}}&=&\dfrac{\partial}{\partial \varepsilon_{t}}(\mu+\varepsilon_{t}+a_{1}\varepsilon_{t-1}+a_{1}^{2}\varepsilon_{t-2}+・・・)\\ &=&1 \end{eqnarray}
となる。①式を1期先にずらすと、

y_{t+1}=\mu+\sum\limits_{i=0}^{\infty} a_{1}^{i}\varepsilon_{t+1-i}
となる。これを\varepsilon_{t}で微分すると、

\begin{eqnarray} \dfrac{\partial y_{t+1}}{\partial \varepsilon_{t}}&=&\dfrac{\partial}{\partial \varepsilon_{t}}(\mu+\varepsilon_{t+1}+a_{1}\varepsilon_{t}+a_{1}^{2}\varepsilon_{t-1}+・・・)\\ &=&a_{1} \end{eqnarray}
となる。①式を2期先にずらすと、

y_{t+2}=\mu+\sum\limits_{i=0}^{\infty} a_{1}^{i}\varepsilon_{t+2-i}
となる。これを\varepsilon_{t}で微分すると、

\begin{eqnarray} \dfrac{\partial y_{t+2}}{\partial \varepsilon_{t}}&=&\dfrac{\partial}{\partial \varepsilon_{t}}(\mu+\varepsilon_{t+2}+a_{1}\varepsilon_{t+1}+a_{1}^{2}\varepsilon_{t}+・・・)\\ &=&a_{1}^{2} \end{eqnarray}
となる。同様に、①式をh期先ずらし、\varepsilon_{t}で微分すると、

\dfrac{\partial y_{t+h}}{\partial \varepsilon_{t}}=a_{1}^{h}・・・②
となる。\dfrac{\partial y_{t+h}}{\partial \varepsilon_{t}}を\theta(h)と書くことにする。この\theta(h)をAR(1)モデルのインパルス応答関数(IRF)とよぶ。

Ⅱ:AR(2)モデルのIRF

AR(p)モデルを考える。(※簡単化のため、定数項は省略する)


y_{t}=a_{1}y_{t-1}+a_{2}y_{t-2}+・・・+a_{p}y_{t-p}+\varepsilon_{t}
AR(p)モデルが定常性を満たしている場合、以下のようなMA(\infty)モデルに書き換えることができる。*2

y_{t}=\sum\limits_{i=0}^{\infty}\theta(i)\varepsilon_{t-i}・・・③
ただし、\theta(i)は\theta(0)=1という初期値から計算を始めて、

\theta(i) = a_{1}\theta(i-1)+a_{2}\theta(i-2)+ ・・・+a_{p}\theta(i-p)
の漸化式に従う係数である。(iが負の場合は\theta(i)=0として計算される。)
p=2として、③式を\varepsilon_{t}で微分すると、

\begin{eqnarray} \dfrac{\partial y_{t}}{\partial \varepsilon_{t}}&=&\dfrac{\partial}{\partial \varepsilon_{t}}(\theta(0)\varepsilon_{t}+\theta(1)\varepsilon_{t-1}+\theta(2)\varepsilon_{t-2}+・・・)\\ &=&\theta(0)\\ &=&1 \end{eqnarray}
となる。②式を1期先にずらして、\varepsilon_{t}で微分すると、

\begin{eqnarray} \dfrac{\partial y_{t+1}}{\partial \varepsilon_{t}}&=&\dfrac{\partial}{\partial \varepsilon_{t}}(\theta(0)\varepsilon_{t+1}+\theta(1)\varepsilon_{t}+\theta(2)\varepsilon_{t-1}+・・・)\\ &=&\theta(1)\\ &=&a_{1}\theta(0)+a_{2}\theta(-1)\\ &=&a_{1} \end{eqnarray}
となる。③式を2期先にずらして、\varepsilon_{t}で微分すると、

\begin{eqnarray} \dfrac{\partial y_{t+2}}{\partial \varepsilon_{t}}&=&\dfrac{\partial}{\partial \varepsilon_{t}}(\theta(0)\varepsilon_{t+2}+\theta(1)\varepsilon_{t+1}+\theta(2)\varepsilon_{t}+・・・)\\ &=&\theta(2)\\ &=&a_{1}\theta(1)+a_{2}\theta(0)\\ &=&a_{1}^{2}+a_{2} \end{eqnarray}
となる。同様に、③式をh期先にずらして、\varepsilon_{t}で微分すると、

\begin{eqnarray} \dfrac{\partial y_{t+h}}{\partial \varepsilon_{t}}&=&\theta(h)\\ &=&a_{1}\theta(h-1)+a_{2}\theta(h-2)・・・④ \end{eqnarray}
となり、④式がAR(2)モデルのIRFである。

Ⅲ:Pythonを用いた実装

Pythonを用いて、AR(1)モデル、AR(2)モデルのIRFを求めてみる。*3以下のライブラリをインポートしておく。

# 必要なライブラリをインポート
import japanize_matplotlib
import numpy as np
import pandas as pd
import statsmodels.api as sm
from matplotlib import pyplot as plt
from matplotlib.pylab import rcParams
rcParams["figure.figsize"] = 15, 6

次にAR(1)モデルのIRFを求めてみる。a_{1}が0.1,0.7,0.95の3パターンで実施する。

#20期、係数として、0.1,0.7,0.95の3パターン
PERIOD = 20
a_1 = [0.1, 0.7, 0.95]

#IRFの定義通りに計算
output = np.zeros([len(a_1), PERIOD])
for i in range(0, len(a_1), 1):
    for j in range(0, PERIOD, 1):
        output[i, j] = 1 * a_1[i] ** j

# グラフの描画
plt.figure(figsize=(10,4))
plt.plot(output[0,:], color='black', label='a_1=0.1')
plt.plot(output[1,:], color='black', linestyle='dashed', label='a_1=0.7')
plt.plot(output[2,:], color='black', linestyle='dotted', label='a_1=0.95')
plt.axhline(0, color="red", linewidth=1)
plt.title('インパルス応答関数', fontsize=18)
plt.xticks(np.arange(0, PERIOD, 1))
plt.xlabel('h', fontsize=14)
plt.ylabel('θ(h)', fontsize=14)
plt.legend(loc='upper right', fontsize=14)

上記のコードを実施した結果が以下である。


図:AR(1)モデルのIRF

a_{1}が0.1,0.7,0.95いずれのパターンでもhが大きくなるにつれて0に近づいていく。しかし、0に収束していく速さは、a_{1}の値が1に近いほど緩やかになる。つまり、a_{1}の値が1に近いほど、ある時点で生じたショックの影響が持続的なものになる。

次に実際のデータを用いて、AR(1)モデル、AR(2)モデルを推定し、推定結果を用いてIRFを求めてみる。使用するのは日本銀行が算出している需給ギャップのデータ*4である。期間は1983年第1四半期～2025年第2四半期である。データをCSVで取得し、データフレームとして取り込み、データを成形する。

#CSVの配置フォルダからCSVを取得
original_data = pd.read_csv(#CSVファイルのパス)
data = original_data.copy()
# "1983.1Q" のような表記を "1983Q1" に整形
data['time'] = data['time'].str.replace('.', 'Q', regex=False)
# "1983.1Q" → "1983Q1" に変換
data['time'] = (
    data['time']
    .astype(str)
    .str.replace(r'\.', 'Q', regex=True)   # 1983.1Q → 1983Q1Q
    .str.replace(r'Q([1-4])Q', r'Q\1', regex=True)  # 1983Q1Q → 1983Q1
    .str.strip()
)
# 四半期末日をDatetimeIndex化
data['date'] = pd.PeriodIndex(data['time'], freq='Q').end_time
data['date'] = data['date'].dt.normalize()  
data = data.drop(columns=['time']).set_index('date')

整形後のデータは以下である。※先頭5行のみ


図:整形後データ

次に1期ラグ、2期ラグ変数を作成する。作成に伴って生じた欠損値が生じるレコードも削除しておく。

# 1期ラグ・2期ラグ列を作成
data["gdpg_lag1"] = data["gdpg"].shift(1)
data["gdpg_lag2"] = data["gdpg"].shift(2)
# 欠損のある最初の2行は削除
data = data.dropna()

ラグ変数の作成と欠損値があるレコードを削除したあとのデータは以下である。※先頭5行のみ


図:ラグ変数作成、欠損値ありレコード削除後データ

以上でデータの準備が完了した。このデータを用いて、AR(1)モデル、AR(2)モデルを推定する。まずはAR(1)モデルを推定する。

# --- AR(1)モデル ---
X1 = sm.add_constant(data["gdpg_lag1"])  # 定数項を追加
y = data["gdpg"]
model_ar1 = sm.OLS(y, X1).fit()
print(model_ar1.summary())

推定結果は以下である。


図:AR(1)モデルの推定結果

次にAR(2)モデルを推定する。

# --- AR(2)モデル ---
X2 = sm.add_constant(data[["gdpg_lag1", "gdpg_lag2"]])  # 定数項を追加
model_ar2 = sm.OLS(y, X2).fit()
print(model_ar2.summary())

推定結果は以下である。


図:AR(2)モデルの推定結果

AR(1)モデル、AR(2)モデルの推定結果を用いて、それぞれのIRFを20期先まで求めていく。与えるショックの大きさを1標準偏差とするため、推定結果から、残差の標準偏差(標準誤差)を取得しておく。

MAX_H = 20  # 20期先まで描画
# --- 標準偏差（残差の標準偏差）を取得 ---
sigma_ar1 = model_ar1.resid.std()
sigma_ar2 = model_ar2.resid.std()

次にIRFを求める。

# --- AR(1)のIRF ---
a_1 = model_ar1.params["gdpg_lag1"]
irf_ar1 = [sigma_ar1]  # t=0 で1標準偏差のショック
for h in range(1, MAX_H):
    irf_ar1.append(irf_ar1[-1] * a_1)#2式に対応している

# --- AR(2)のIRF ---
a_1, a_2 = model_ar2.params["gdpg_lag1"], model_ar2.params["gdpg_lag2"]
irf_ar2 = [sigma_ar2, a_1 * sigma_ar2]  # t=0 で1標準偏差のショック、a_1 * sigma_ar2はt=1のときのIRFの値
for h in range(2, MAX_H):
    irf_ar2.append(a_1 * irf_ar2[-1] + a_2 * irf_ar2[-2])#4式に対応している

# --- グラフの描画 ---
plt.figure(figsize=(10, 5))
plt.plot(irf_ar1, label="AR(1)", linewidth=2, color="black")
plt.plot(irf_ar2, label="AR(2)", linewidth=2, linestyle="dashed", color="gray")
plt.axhline(0, color="red", linewidth=1)
plt.xticks(np.arange(0, len(irf_ar1), 1))
plt.title("インパルス応答関数（1標準偏差ショック）", fontsize=16)
plt.xlabel("四半期", fontsize=12)
plt.legend(fontsize=12)
plt.grid(False)
plt.show()

求めたIRFを可視化した結果が以下である。


図:AR(1)モデル、AR(2)モデルのIRF

AR(1)モデルのIRFは単調に減少している。一方、AR(2)モデルのIRFはショックが発生した数期後に最大値を取り、そこから減少している。このような形状のIRFを「こぶ型」とよぶ。最後に、AICを算出し、AR(1)モデルとAR(2)モデルのどちらが適しているかを判定する。

# AIC値の取得
aic_ar1 = model_ar1.aic
aic_ar2 = model_ar2.aic

# 比較出力
print("=== AICによるモデル比較 ===")
print(f"AR(1)モデルのAIC: {aic_ar1:.3f}")
print(f"AR(2)モデルのAIC: {aic_ar2:.3f}")

if aic_ar1 < aic_ar2:
    print("→ AICに基づいてラグ次数を選択するとAR(1)が選択される")
else:
    print("→ AICに基づいてラグ次数を選択するとAR(2)が選択される")

両モデルについてAICを算出した結果が以下である。


図:両モデルのAIC

AICの結果からはAR(2)モデルが望ましいという結果が得られた。IRFもAR(2)モデルのものが望ましい。

Ⅳ:参考文献

●雑誌
新谷元嗣・前橋昂平(2024)「マクロ経済政策評価のための時系列分析　第1回」『経済セミナー』2024年10・11月号,p49-58,日本評論社
●Webサイト
経済セミナー連載「マクロ経済政策評価のための時系列分析」サポートサイト

*1:過去記事を参照してください。

*2:新谷・前橋[2024]を参照してください。ここの導出は稿を改めて書きたいと思っています。

*3:以下のPythonコードは新谷・前橋[2024]のサポートサイトで公開されているものをベースにして、本稿用に変更を加えたものです。

*4:ここから取得しました。

身内用

ライブラリインポート

#基本ライブラリ
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
#ダウンサンプリング
from imblearn.under_sampling import RandomUnderSampler, ClusterCentroids, NearMiss
#標準化、ワンホットエンコーディング、順序化
from sklearn.preprocessing import StandardScaler, OneHotEncoder, OrdinalEncoder
#クラスタリング、決定木
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
from sklearn.tree import DecisionTreeClassifier, plot_tree
from sklearn.model_selection import train_test_split
from sklearn.metrics import classification_report
#主成分分析
from sklearn.decomposition import PCA
import seaborn as sns
#ロジスティック回帰
import statsmodels.api as sm
#その他
import japanize_matplotlib
import os
import warnings
warnings.simplefilter('ignore')

サンプルデータ生成

np.random.seed(42)  
# サンプルサイズ
n = 1000
data = {
    "顧客ID": np.arange(1, n+1),
    "苦情フラグ": np.random.choice([0, 1], size=n, p=[0.8, 0.2]),
    "年代": np.random.choice(["10代", "20代", "30代", "40代"], size=n),
    "性別": np.random.choice(["男性", "女性"], size=n),
    "居住地域": np.random.choice(["関東", "関西", "東北", "その他"], size=n),
    "契約回数": np.random.randint(1, 16, size=n),
    "職員フラグ": np.random.choice([0, 1], size=n), 
    "職員役職": np.random.choice(["新人", "中堅", "ベテラン"], size=n),
}

imported_df = pd.DataFrame(data).set_index("顧客ID")

データ読み込み

# カレントディレクトリを基準に "分析データ" フォルダを指定
data_dir = os.path.join(os.getcwd(), "分析データ")
# CSVファイル名を指定
file_name = "data.csv"
# フルパスを生成
file_path = os.path.join(data_dir, file_name)

# CSVの読み込み処理
if mode == "all":
    imported_df = pd.read_csv(file_path)  # 全列を取り込む
elif mode == "subset":
    imported_df = pd.read_csv(file_path, usecols=selected_columns)  # 指定列だけを取り込む
else:
    raise ValueError("mode は 'all' か 'subset' を指定してください")

#顧客番号をindex化する
imported_df = pd.DataFrame(data).set_index("顧客番号")  

#コピー
df = imported_df.copy()

データ前処理

# ==========================
# 前処理を実施する列を指定する
# ==========================
# ①そのまま
cols_keep = [
    "苦情フラグ",
    "退社フラグ"
]  
# ②標準化対象カラム
cols_standardize = [
    "契約回数"
]
# ③順序ありの数値変換対象カラム
cols_ordinal = [
    "役職", 
    "等級"
] 
# ④順序なしカテゴリ化対象カラム
cols_nominal = [
    "年代", 
    "性別", 
    "居住地域"
]  

# ==========================
# 前処理を実施
# ==========================
#標準化
scaler = StandardScaler()
df_standardized = pd.DataFrame(
    scaler.fit_transform(df[cols_standardize]),
    columns=cols_standardize,
    index=df.index
)
#順序の数値化
ordinal_mappings = [
    #※cols_ordinalの記載順と対応させる。
    ["rankA", "rnakB", "rankC"],  
    ["rank1", "rank2", "rank3"]   

] 
ord_encoder = OrdinalEncoder(categories=ordinal_mappings)
df_ordinal = pd.DataFrame(
    ord_encoder.fit_transform(df[cols_ordinal]),
    columns=cols_ordinal,
    index=df.index
)
#順序なしカテゴリ化
nominal_encoder = OneHotEncoder(drop=None, sparse_output=True)#※クラスタリングのときはNone、ロジ回のときはfirst
nominal_sparse = nominal_encoder.fit_transform(df[cols_nominal])
df_nominal_ohe = pd.DataFrame(
    nominal_sparse.toarray(),
    columns=nominal_encoder.get_feature_names_out(cols_nominal),
    index=df.index
)
#前処理不要列
df_keep = df[cols_keep] 
#すべて結合
df_processed = pd.concat(
    [df_keep, df_standardized, df_ordinal, df_nominal_ohe],
    axis=1
)

ダウンサンプリング

# ==========================
# 苦情フラグカラムとそれ以外のカラムを分割
# ==========================
Y=df_processed["苦情フラグ"]
X=df_processed.drop(columns=["苦情フラグ"])

# ==========================
# 1. RandomUnderSampler
# ==========================
rus = RandomUnderSampler(sampling_strategy=1.0, random_state=42)
X_rus, y_rus = rus.fit_resample(X, Y)

df_rus = pd.DataFrame(X_rus, columns=X.columns)
df_rus["苦情フラグ"] = y_rus

print("\n=== RandomUnderSampler 後 ===")
print(df_rus["苦情フラグ"].value_counts().to_dict())

# ==========================
# 2. ClusterCentroids
# ==========================
cc = ClusterCentroids(sampling_strategy=1.0, random_state=42)
X_cc, y_cc = cc.fit_resample(X, Y)

df_cc = pd.DataFrame(X_cc, columns=X.columns)
df_cc["苦情フラグ"] = y_cc

print("\n=== ClusterCentroids 後 ===")
print(df_cc["苦情フラグ"].value_counts().to_dict())

# ==========================
# 3. NearMiss
# ==========================
nm = NearMiss(version=1)
X_nm, y_nm = nm.fit_resample(X, Y)
df_nm = pd.DataFrame(X_nm, columns=X.columns)
df_nm["苦情フラグ"] = y_nm

print("\n=== NearMiss 後 ===")
print(df_nm["苦情フラグ"].value_counts().to_dict())


クラスタリング

# 特徴量と目的変数
X_rus = df_rus.drop(columns=["苦情フラグ"])
y_rus = df_rus["苦情フラグ"]

# ==========================
# k-meansクラスタリング実施関数
# ==========================
def kmeans_clustering(
    df,bef_df, label_col="苦情フラグ", K=range(2,6), random_state=42
):
    inertia, silhouettes = [], []
    K_list = list(K)
    X = df.drop(columns=["苦情フラグ"])

    # ==========================
    # 1. エルボー & シルエット
    # ==========================
    for k in K_list:
        km = KMeans(n_clusters=k, random_state=random_state, n_init="auto")
        labels = km.fit_predict(X)
        inertia.append(km.inertia_)
        silhouettes.append(silhouette_score(X, labels))

    # 図
    fig, ax = plt.subplots(1, 2, figsize=(14, 5))
    ax[0].plot(K_list, inertia, marker="o")
    ax[0].set_title("エルボー法によるクラスタ数の決定")
    ax[0].set_xlabel("クラスタ数 (k)")
    ax[0].set_ylabel("SSE (Inertia)")
    ax[0].grid(True)

    ax[1].plot(K_list, silhouettes, marker="o")
    ax[1].set_title("シルエットスコアによるクラスタ数の評価")
    ax[1].set_xlabel("クラスタ数 (k)")
    ax[1].set_ylabel("シルエットスコア")
    ax[1].grid(True)
    plt.show()

    # ==========================
    # 2. 最適kを選択
    # ==========================
    best_k = K_list[int(np.argmax(silhouettes))]
    print("シルエットスコア最大のクラスタ数:", best_k)

    # ==========================
    # 3. KMeans再実行 & クラスタ付与
    # ==========================
    km_final = KMeans(n_clusters=best_k, random_state=random_state, n_init="auto")
    clusters = km_final.fit_predict(X)
    df_result = df.copy()
    df_result["クラスタ"] = clusters

    # ==========================
    # 4. 件数・苦情率
    # ==========================
    cluster_size = df_result["クラスタ"].value_counts().sort_index()
    cluster_label_mean = df_result.groupby("クラスタ")[label_col].mean()

    print("\n=== クラスタごとの件数 ===")
    print(cluster_size)
    print("\n=== クラスタごとの{}平均 ===".format(label_col))
    print(cluster_label_mean)
    return df_result, cluster_size, cluster_label_mean

決定木

# 決定木を用いてクラスタ分類ルールを学習・可視化する関数
# Parameters
# ----------
# df_rus : pd.DataFrame
#     「クラスタ」と「苦情フラグ」を含むデータフレーム
# test_size : float, default=0.3
#     train/test 分割のテストデータ比率
# random_state : int, default=42
#     乱数シード
# max_depth : int, default=4
#     決定木の最大深さ
# Returns
# -------
# dt : DecisionTreeClassifier
#     学習済み決定木モデル
# feature_importance : pd.Series
#     特徴量重要度
def decision_tree_cluster( df_rus,test_size=0.3,random_state=42,max_depth=4):
    # ==========================================
    # 説明変数（苦情フラグとクラスタは除外）
    # ==========================================
    X_tree = df_rus.drop(columns=["クラスタ", "苦情フラグ"])
    y_tree = df_rus["クラスタ"]

    # train/test 分割
    X_train, X_test, y_train, y_test = train_test_split(
        X_tree, y_tree,
        test_size=test_size,
        random_state=random_state,
        stratify=y_tree
    )

    # ==========================================
    # 決定木モデルの学習
    # ==========================================
    dt = DecisionTreeClassifier(max_depth=max_depth, random_state=random_state)
    dt.fit(X_train, y_train)

    # ==========================================
    # 精度評価
    # ==========================================
    y_pred = dt.predict(X_test)
    print("=== 決定木によるクラスタ分類 ===")
    print(classification_report(y_test, y_pred))

    # ==========================================
    # 特徴量重要度
    # ==========================================
    feature_importance = pd.Series(dt.feature_importances_, index=X_tree.columns)
    print("\n=== 特徴量重要度 ===")
    print(feature_importance.sort_values(ascending=False))

    feature_importance.sort_values(ascending=True).plot(kind="barh", figsize=(8,5))
    plt.title("決定木に基づく特徴量重要度")
    plt.show()

    # ==========================================
    # 決定木の可視化
    # ==========================================
    plt.figure(figsize=(14,8))
    plot_tree(
        dt,
        feature_names=X_tree.columns,
        class_names=[str(c) for c in sorted(y_tree.unique())],
        filled=True, rounded=True, fontsize=10
    )
    plt.title("クラスタ分類を説明する決定木")
    plt.show()
    return dt, feature_importance

ロジスティック回帰

y_logi = df_processed_logi["苦情フラグ"]
X_logi = df_processed_logi.drop(columns=["苦情フラグ"])

# 定数項を追加
X = sm.add_constant(X)

# ===== 3. ロジスティック回帰の実行 =====
model = sm.Logit(y_logi, X_logi)
result = model.fit()

# ===== 4. 結果の表示 =====
print(result.summary())

# ===== 5. オッズ比表を作成 =====
odds_ratios = pd.DataFrame({
    "変数": result.params.index,
    "係数": result.params.values,
    "オッズ比": np.exp(result.params.values),
    "p値": result.pvalues.values
})
print(odds_ratios)

経済学メモ:経済学でよく使う数式「等差・等比数列と単利・複利、割引現在価値」

経済学 経済学でよく使う数式

◆本稿の内容
等差・等比数列の一般項(等比数列については和も)を確認した後に、それぞれの応用例として、資金の単利運用、複利運用による元利合計額の計算方法を紹介します。投資の世界でよく知られている「72ルール」にも触れています。等比数列、複利計算の更なる応用として、割引現在価値の概念も紹介します。

◆目次

Ⅰ.等差数列と等比数列

Ⅰ-1.等差数列

初項a、公差dの等差数列\{a_{t}\}は次のように書ける。


a_{0}=a\\ a_{1}=a+d\\ a_{2}=a+2d\\ \vdots\\ a_{t}=a+td\\ \vdots

よって、初項a、公差dの等差数列\{a_{t}\}の一般項は


a_{t}=a+td　(t=0,1,2,・・・)

となる。*1

Ⅰ-2.等比数列

初項a、公比rの等比数列\{a_{t}\}は次のように書ける。


a_{0}=a\\ a_{1}=ar\\ a_{2}=ar^{2}\\ \vdots\\ a_{t}=ar^{t}\\ \vdots

よって、初項a、公比rの等比数列\{a_{t}\}の一般項は


a_{t}=ar^{t}　(t=0,1,2,・・・)

となる。*2等比数列の初項a(0番目の値)からt番目までの項の和S_{t}の公式を求めておく。S_{t}は


\begin{eqnarray} S_{t}&=&a+a_{1}+a_{2}+\dots +a_{t}\\ &=&a+ar+ar^{2}+\dots +ar^{t} \end{eqnarray}

となる。この式からrS_{t}=ar+ar^{2}++ar^{3+}\dots +ar^{t}+ar^{t+1}を引いて式を整理すると、


\begin{eqnarray} S_{t}-rS_{t}&=&a-ar^{t+1}\\ (1-r)S_{t}&=&a(1-r^{t+1})\\ S_{t}&=&\dfrac{a(1-r^{t+1})}{1-r} \end{eqnarray}

となる。S_{t}のtを\inftyとして、無限和S_{\infty}を考える。|r|が1未満ならば、r^{t+1}の極限は0となる。よって、S_{\infty}は\dfrac{a}{1-r}となる。S_{\infty}は「無限等比級数の和」と呼ぶ。

Ⅱ.単利と等差数列、複利と等比数列

Ⅱ-1:単利と等差数列

単利とは、初期時点の元金のみを対象にして利息を計算する方法である。例として、元金10000円を年利5%で、3年間単利で運用する場合の元利合計額を考える。まず、毎年ごとに得られる利息を次のように求める。


10000×0.05=500

3年間運用するから、1年間分の利息(500円)を3回得られることになる。元金に3年間の運用で得た利息の合計を足し合わせることで元利合計を次のように求めることができる。


10000+500×3=11500

この式と等差数列の一般項との対応関係を考えてみる。元金の10000円は等差数列の一般項の初項aに、1年間分の利息の500円は公差dに対応している。つまり、等差数列の一般項の形で表すと、


a_{t}=10000+500t　(t=0,1,2,・・・)

と書ける。上式にt=3を代入すると3年後の元利合計を求めることができる。

Ⅱ-2.複利と等比数列

複利とは、前期時点までの元利合計を対象にして今期の利息を計算する方法である。例として、元金10000円を年利5%で、3年間複利で運用する場合の元利合計額を考える。1年後の元利合計は


10000×(1+0.05)=10500

となる。2年後の元利合計はこの10500円に対して5%の利息が付く形になるから、


10500×(1+0.05)=11025

となる。3年後の元利合計はこの11025円に対して5%の利息が付く形になるから、


11025×(1+0.05)≒11576

となる。この計算過程を前提として上式を


10000×(1+0.05)^{3}≒11576

と書き換えることができる。この式と等比数列の一般項との対応関係を考えてみる。元金の10000円は等比数列の一般項の初項aに、1に年利を足した値(1+0.05)*3が公比rに対応している。つまり、等比数列の一般項の形で表すと


a_{t}=10000×(1+0.05)^{t}　(t=0,1,2,・・・)

と書ける。上式にt=3を代入すると3年後の元利合計を求めることができる。

Ⅱ-3.元金が2倍になるまでの期間:72ルール

複利運用の場合に元金が2倍になるまでの期間を求めてみよう。元金をa、年利をr\%とする。等比数列の一般項の左辺に元金aを2倍した値を代入して式を整理する。


\begin{eqnarray} a_{t}&=&a\left(1+\dfrac{r}{100}\right)^{t}　(t=0,1,2,・・・)\\ 2a&=&a\left(1+\dfrac{r}{100}\right)^{t}\\ 2&=&\left(1+\dfrac{r}{100}\right)^{t} \end{eqnarray}

両辺の自然対数を取り、tについて解くと、


\begin{eqnarray} \log 2 &=& \log \left(1+\dfrac{r}{100}\right)^{t}\\ \log 2 &=&t \log \left(1+\dfrac{r}{100}\right)\\ t&=&\dfrac{\log 2}{\log \left(1+\dfrac{r}{100}\right)} \end{eqnarray}

となる。ここでxが小さいとき、\log (1+x)はxで近似できることを利用すると、*4


\begin{eqnarray} t&=&\dfrac{\log 2}{\log \left(1+\dfrac{r}{100}\right)}\\ t&=&\dfrac{\log 2}{\dfrac{r}{100}}\\ \end{eqnarray}

となる。この式を用いて、年利がr\%のときに元金が何年で2倍になるかを計算した結果が以下の表である。


図1:2倍になるまでの年数

1列目が年利、2列目が2倍になるまでの年数である。3列目は72でrを割った値である。2列目と3列目を比較すると、おおよそ等しい値になっていることが分かる。\log 2はおよそ0.6931である。だが、暗算のときに計算がしづらい。そこでに0.6931に近く、約数が多く暗算のしやすい0.72を用いて、\dfrac{0.72}{\dfrac{r}{100}}=\dfrac{72}{r}という式で、元金が2倍になるまでのおおよその年数を計算することがある。これを「72ルール」や「72の法則」と呼ぶ。

Ⅲ:割引現在価値

Ⅲ-1.現在価値に割り引くとは

年利r\%で資金を運用する手段があるとする。このとき、t年後に得られる金額の価値を現在時点で評価するとどのような値になるだろうか。c円をこの運用手段でt年運用するとt年後の元利合計はc\left(1+\dfrac{r}{100}\right)^{t}になる。この値をaとおいてcについて解くと、


\begin{eqnarray} c\left(1+\dfrac{r}{100}\right)^{t}&=&a\\ c&=&\dfrac{a}{\left(1+\dfrac{r}{100}\right)^{t}}\\ \end{eqnarray}

となる。この式はt年後に得られるa円は現在時点のc=\dfrac{a}{\left(1+\dfrac{r}{100}\right)^{t}}円に等しいということを示している。この値をt年後に得られるa円の割引現在価値と呼ぶ。このように、将来時点で得られる金額を現在時点の金額(価値)に割り引いて考えることを「現在価値に割り引く」と呼ぶ。割引の際に用いる因子(上式でいうと、1+\dfrac{r}{100}のこと)を割引因子と呼ぶ。

Ⅲ-2.割引現在価値の総和

割引現在価値の応用例として、投資プロジェクトの評価を考える。今、初期時点において1000万円を投資すると、1年後から毎年40万円の収益を継続して得られる*5プロジェクトがあるとする。ただし、割引因子として使用する金利を3%とする。このプロジェクトは費用対効果の観点で実施すべきであろうか。この問題を考えるにあたって、まずはこの投資プロジェクトの費用を考える。費用は初期時点においてのみ必要な1000万円である。次に収益を考える。1年後から継続して毎年40万円の収益を得られるから、この毎年得られる収益を現在価値に割り引いて考える必要がある。1年後に得られる収益を割引現在価値に直すと\dfrac{40}{1+0.03}円である。2年後に得られる収益を割引現在価値に直すと\dfrac{40}{(1+0.03)^{2}}円である。3年後に得られる収益を割引現在価値に直すと\dfrac{40}{(1+0.03)^{3}}円である。同様に考えていくと、t年後に得られる収益の割引現在価値は\dfrac{40}{(1+0.03)^{t}}円である。これらの割引現在価値の総和がこの投資プロジェクトの割引現在価値である。割引現在価値の総和PVは


\begin{eqnarray} PV&=&\dfrac{40}{1+0.03}+\dfrac{40}{(1+0.03)^{2}}+\dfrac{40}{(1+0.03)^{3}}+\dots+\dfrac{40}{(1+0.03)^{t}}+\dots\\ &=&\sum_{t=1}^{\infty}\dfrac{40}{(1+0.03)^{t}}\\ &=&\dfrac{\dfrac{40}{1+0.03}}{1-\dfrac{1}{1+0.03}}\\ &=&\dfrac{40}{0.03}\\ &≒&1333 \end{eqnarray}

となる。*6よって、
「プロジェクト実行により生じる現在時点で必要な費用」<「プロジェクト実行によって将来に渡って得られる収益の割引現在価値の総和」
となるから、この投資プロジェクトは実行すべきである。*7

Ⅳ:参考文献

尾山大輔・安田洋祐編著(2013)『改訂版 経済学で出る数学ー高校数学からきちんと攻める』日本評論社
高木悟・上江洲弘明(2022)『早稲田大学全学教育基盤シリーズ 金利の計算 解析学への入り口』共立出版

*1:高校数学では初項の番号を1にしていますが、本稿では0にしています。この違いにより、等差数列の一般項の形も高校数学で出てくるa_{t}=a+(t-1)d　(t=1,2,・・・)とは少し異なっています。この違いは見た目上のものであって、本質的な違いは全くありません。

*2:ここでも初項の番号を0にしています。

*3:利子率に1を足した値のことを「粗」利子率と呼びます。利子率は粗利子率との対比を意識して「純」利子率とも呼びます。

*4:これを参照してください。

*5:初期時点の投資のみで永遠に収益を上げ続けることのできるプロジェクトは存在しないので現実離れしている設定だと思いますが気にしないでください。無限期間における割引現在価値の総和の計算方法を確認することが主目的です。

*6:2行目から3行目の式変形で無限等比級数の和の式を使っています。

*7:4行目の式から、投資プロジェクトの毎期の収益が大きくなるほど、割引因子として使用する金利が小さくなるほど、投資プロジェクトの割引現在価値の総和が大きくなり、実行されやすくなることが分かります。