Progress. The TODOs still remain

Author: gusamarante

example_us.py

@@ -16,5 +16,28 @@
     selected_maturities=[6, 12, 24, 36, 48, 60, 72, 84, 96, 108, 120]
 )
 
-acm.tp[119].plot()
+# CHART
+size = 7
+fig = plt.figure(figsize=(size * (16 / 7.3), size))
+
+ax = plt.subplot2grid((1, 2), (0, 0))
+ax.plot(ylds_d[120], label="Actual Yield", lw=1)
+ax.plot(acm.miy[120], label="Fitted Yield", lw=1, ls='--')
+ax.set_title("10-Year Model Fit")
+ax.xaxis.grid(color="grey", linestyle="-", linewidth=0.5, alpha=0.5)
+ax.yaxis.grid(color="grey", linestyle="-", linewidth=0.5, alpha=0.5)
+ax.tick_params(rotation=90, axis="x")
+ax.legend(loc="upper right")
+
+ax = plt.subplot2grid((1, 2), (0, 1))
+ax.plot(ylds_d[120], label="Yield", lw=1)
+ax.plot(acm.rny[120], label="Risk Neutral Yield", lw=1)
+ax.plot(acm.tp[120], label="Term Premium", lw=1)
+ax.set_title("10-Year Yield Decomposition")
+ax.xaxis.grid(color="grey", linestyle="-", linewidth=0.5, alpha=0.5)
+ax.yaxis.grid(color="grey", linestyle="-", linewidth=0.5, alpha=0.5)
+ax.tick_params(rotation=90, axis="x")
+ax.legend(loc="upper right")
+
+plt.tight_layout()
 plt.show()

pyacm/acm.py

@@ -3,11 +3,13 @@
 
 from numpy.linalg import inv
 from sklearn.decomposition import PCA
+from statsmodels.tools.tools import add_constant
 
 from pyacm.utils import vec, vec_quad_form, commutation_matrix
 
 # TODO Curve in daily frequency could be None
 # TODO Make sure it works for DI Futures
+# TODO make sure data is accepted as decimals
 class NominalACM:
     """
     This class implements the model from the article:
@@ -137,6 +139,7 @@ def __init__(
 
         # TODO assert columns of daily and monthly are the same
         # TODO assert monthly index frequency
+        # TODO assert columns are consecutive integers
 
         self.n_factors = n_factors
         self.curve = curve
@@ -163,21 +166,36 @@ def __init__(
 
         # 3rd Step - Convexity-adjusted price of risk
         self.lambda0, self.lambda1, self.mu_star, self.phi_star = self._retrieve_lambda()
-        # TODO EVERYTHING RIGHT UP TO HERE
 
-        # TODO PAREI AQUI - Olhar em que situações pode-se usar os yield mensais.
-        if self.curve.index.freqstr == 'M':
-            X = self.pc_factors_m
-            r1 = self.curve_monthly.iloc[:, 0]
-        else:
-            X = self.pc_factors_d
-            r1 = self.curve.iloc[:, 0]
+        # Short Rate Equation
+        self.delta0, self.delta1 = self._short_rate_equation(
+            r1=self.curve_monthly.iloc[:, 0],
+            X=self.pc_factors_m,
+        )
+
+        # Affine Yield Coefficients
+        self.A, self.B = self._affine_coefficients(
+            lambda0=self.lambda0,
+            lambda1=self.lambda1,
+        )
+
+        # Risk-Neutral Coefficients
+        self.Arn, self.Brn = self._affine_coefficients(
+            lambda0=np.zeros(self.lambda0.shape),
+            lambda1=np.zeros(self.lambda1.shape),
+        )
 
-        self.miy = self._affine_recursions(self.lambda0, self.lambda1, X, r1)
-        self.rny = self._affine_recursions(0, 0, X, r1)
+        # Model Implied Yield
+        self.miy = self._compute_yields(self.A, self.B)
+
+        # Risk Neutral Yield
+        self.rny = self._compute_yields(self.Arn, self.Brn)
+
+        # Term Premium
         self.tp = self.miy - self.rny
-        self.er_loadings, self.er_hist_m, self.er_hist_d = self._expected_return()
-        self.z_lambda, self.z_beta = self._inference()
+
+        # Expected Return
+        self.er_loadings, self.er_hist = self._expected_return()
 
     def fwd_curve(self, date=None):
         """
@@ -340,34 +358,43 @@ def _retrieve_lambda(self):
 
         return lambda0, lambda1, muStar, phiStar
 
-    def _affine_recursions(self, lambda0, lambda1, X_in, r1):
-        X = X_in.T.values[:, 1:]
-        r1 = vec(r1.values)[-X.shape[1]:, :]
-
-        A = np.zeros((1, self.n))
-        B = np.zeros((self.n_factors, self.n))
-
-        delta = r1.T @ np.linalg.pinv(np.vstack((np.ones((1, X.shape[1])), X)))
-        delta0 = delta[[0], [0]]
-        delta1 = delta[[0], 1:]
-
-        A[0, 0] = - delta0
-        B[:, 0] = - delta1
-
-        for i in range(self.n - 1):
-            A[0, i + 1] = A[0, i] + B[:, i].T @ (self.mu - lambda0) + 1 / 2 * (B[:, i].T @ self.Sigma @ B[:, i] + 0 * self.sigma2) - delta0
-            B[:, i + 1] = B[:, i] @ (self.phi - lambda1) - delta1
-
-        # Construct fitted yields
-        ttm = np.arange(1, self.n + 1) / 12
-        fitted_log_prices = (A.T + B.T @ X).T
-        fitted_yields = - fitted_log_prices / ttm
-        fitted_yields = pd.DataFrame(
-            data=fitted_yields,
-            index=self.curve.index[1:],
+    @staticmethod
+    def _short_rate_equation(r1, X):
+        r1 = r1 / 1200  # TODO remove the 100
+        X = add_constant(X)
+        Delta = inv(X.T @ X) @ X.T @ r1
+        delta0 = Delta.iloc[0]
+        delta1 = Delta.iloc[1:].values
+        return delta0, delta1
+
+    def _affine_coefficients(self, lambda0, lambda1):
+        lambda0 = lambda0.reshape(-1, 1)
+
+        A = np.zeros(self.n)
+        B = np.zeros((self.n, self.n_factors))
+
+        A[0] = - self.delta0
+        B[0, :] = - self.delta1
+
+        for n in range(1, self.n):
+            Bpb = np.kron(B[n - 1, :], B[n - 1, :])
+            s0term = 0.5 * (Bpb @ self.s0 + self.omega[0, 0])
+
+            A[n] = A[n - 1] + B[n - 1, :] @ (self.mu - lambda0) + s0term + A[0]
+            B[n, :] = B[n - 1, :] @ (self.phi - lambda1) + B[0, :]
+
+        return A, B
+
+    def _compute_yields(self, A, B):
+        A = A.reshape(-1, 1)
+        multiplier = np.tile(self.curve.columns / 12, (self.t_d, 1)).T / 100  # TODO remove the 100
+        yields = (- ((np.tile(A, (1, self.t_d)) + B @ self.pc_factors_d.T) / multiplier).T).values
+        yields = pd.DataFrame(
+            data=yields,
+            index=self.curve.index,
             columns=self.curve.columns,
         )
-        return fitted_yields
+        return yields
 
     def _expected_return(self):
         """
@@ -378,93 +405,20 @@ def _expected_return(self):
         expected returns
         """
         stds = self.pc_factors_m.std().values[:, None].T
-        er_loadings = (self.beta.T @ self.lambda1) * stds
+        er_loadings = (self.B @ self.lambda1) * stds
         er_loadings = pd.DataFrame(
             data=er_loadings,
             columns=self.pc_factors_m.columns,
-            index=self.selected_maturities,
+            index=range(1, self.n + 1),
         )
 
-        # Monthly
-        exp_ret = (self.beta.T @ (self.lambda1 @ self.pc_factors_m.T + self.lambda0)).values
-        conv_adj = np.diag(self.beta.T @ self.Sigma @ self.beta) + self.sigma2
-        er_hist = (exp_ret + conv_adj[:, None]).T
-        er_hist_m = pd.DataFrame(
-            data=er_hist,
-            index=self.pc_factors_m.index,
-            columns=self.curve.columns[:er_hist.shape[1]]
-        )
-
-        # Higher frequency
-        exp_ret = (self.beta.T @ (self.lambda1 @ self.pc_factors_d.T + self.lambda0)).values
-        conv_adj = np.diag(self.beta.T @ self.Sigma @ self.beta) + self.sigma2
+        # Historical estimate
+        exp_ret = (self.B @ (self.lambda1 @ self.pc_factors_d.T + self.lambda0.reshape(-1, 1))).values
+        conv_adj = np.diag(self.B @ self.Sigma @ self.B.T) + self.omega[0, 0]
         er_hist = (exp_ret + conv_adj[:, None]).T
         er_hist_d = pd.DataFrame(
             data=er_hist,
             index=self.pc_factors_d.index,
-            columns=self.curve.columns[:er_hist.shape[1]]
+            columns=self.curve.columns,
         )
-
-        return er_loadings, er_hist_m, er_hist_d
-
-    def _inference(self):
-        # TODO I AM NOT SURE THAT THIS SECTION IS CORRECT
-
-        # Auxiliary matrices
-        Z = self.pc_factors_m.copy().T
-        Z = Z.values[:, 1:]
-        Z = np.vstack((np.ones((1, self.t)), Z))
-
-        Lamb = np.hstack((self.lambda0, self.lambda1))
-
-        rho1 = np.zeros((self.n_factors + 1, 1))
-        rho1[0, 0] = 1
-
-        A_beta = np.zeros((self.n_factors * self.beta.shape[1], self.beta.shape[1]))
-
-        for ii in range(self.beta.shape[1]):
-            A_beta[ii * self.beta.shape[0]:(ii + 1) * self.beta.shape[0], ii] = self.beta[:, ii]
-
-        BStar = np.squeeze(np.apply_along_axis(vec_quad_form, 1, self.beta.T))
-
-        comm_kk = commutation_matrix(shape=(self.n_factors, self.n_factors))
-        comm_kn = commutation_matrix(shape=(self.n_factors, self.beta.shape[1]))
-
-        # Assymptotic variance of the betas
-        v_beta = self.sigma2 * np.kron(np.eye(self.beta.shape[1]), inv(self.Sigma))
-
-        # Assymptotic variance of the lambdas
-        upsilon_zz = (1 / self.t) * Z @ Z.T
-        v1 = np.kron(inv(upsilon_zz), self.Sigma)
-        v2 = self.sigma2 * np.kron(inv(upsilon_zz), inv(self.beta @ self.beta.T))
-        v3 = self.sigma2 * np.kron(Lamb.T @ self.Sigma @ Lamb, inv(self.beta @ self.beta.T))
-
-        v4_sim = inv(self.beta @ self.beta.T) @ self.beta @ A_beta.T
-        v4_mid = np.kron(np.eye(self.beta.shape[1]), self.Sigma)
-        v4 = self.sigma2 * np.kron(rho1 @ rho1.T, v4_sim @ v4_mid @ v4_sim.T)
-
-        v5_sim = inv(self.beta @ self.beta.T) @ self.beta @ BStar
-        v5_mid = (np.eye(self.n_factors ** 2) + comm_kk) @ np.kron(self.Sigma, self.Sigma)
-        v5 = 0.25 * np.kron(rho1 @ rho1.T, v5_sim @ v5_mid @ v5_sim.T)
-
-        v6_sim = inv(self.beta @ self.beta.T) @ self.beta @ np.ones((self.beta.shape[1], 1))
-        v6 = 0.5 * (self.sigma2 ** 2) * np.kron(rho1 @ rho1.T, v6_sim @ v6_sim.T)
-
-        v_lambda_tau = v1 + v2 + v3 + v4 + v5 + v6
-
-        c_lambda_tau_1 = np.kron(Lamb.T, inv(self.beta @ self.beta.T) @ self.beta)
-        c_lambda_tau_2 = np.kron(rho1, inv(self.beta @ self.beta.T) @ self.beta @ A_beta.T @ np.kron(np.eye(self.beta.shape[1]), self.Sigma))
-        c_lambda_tau = - c_lambda_tau_1 @ comm_kn @ v_beta @ c_lambda_tau_2.T
-
-        v_lambda = v_lambda_tau + c_lambda_tau + c_lambda_tau.T
-
-        # extract the z-tests
-        sd_lambda = np.sqrt(np.diag(v_lambda).reshape(Lamb.shape, order='F'))
-        sd_beta = np.sqrt(np.diag(v_beta).reshape(self.beta.shape, order='F'))
-
-        z_beta = pd.DataFrame(self.beta / sd_beta, index=self.pc_factors_m.columns, columns=self.selected_maturities).T
-        z_lambda = pd.DataFrame(Lamb / sd_lambda, index=self.pc_factors_m.columns, columns=[f"lambda {i}" for i in range(Lamb.shape[1])])
-
-        return z_lambda, z_beta
-
-
+        return er_loadings, er_hist_d

setup.py

@@ -26,6 +26,7 @@
         'numpy',
         'pandas',
         'scikit-learn',
+        'statsmodels',
         'tqdm',
     ],
     keywords=[