From e1645b7c507ce188501a28a49f2c5452b3e468a5 Mon Sep 17 00:00:00 2001 From: lianghx123 <68636898+lianghx123@users.noreply.github.com> Date: Sun, 30 Jun 2024 01:38:18 +0800 Subject: [PATCH 1/5] grid search of starting_values is not correct --- src/frds/algorithms/_mgarch.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/frds/algorithms/_mgarch.py b/src/frds/algorithms/_mgarch.py index c2c8a14..f38fdf3 100644 --- a/src/frds/algorithms/_mgarch.py +++ b/src/frds/algorithms/_mgarch.py @@ -313,6 +313,7 @@ def starting_values(self) -> Tuple[float, float]: continue ll = -self.loglikelihood_model([a, b]) if ll > max_ll: + max_ll = ll initial_values = [a, b] return initial_values From b70724e88ba66b5e81e4b800fb3c59beaf25045b Mon Sep 17 00:00:00 2001 From: lianghx123 <68636898+lianghx123@users.noreply.github.com> Date: Sun, 30 Jun 2024 01:52:09 +0800 Subject: [PATCH 2/5] Fix missing fit state check in GJRGARCHModel.fit() Added a crucial state check to the fit method in GJRGARCHModel that was previously omitted. The absence of this check led to redundant estimations during the optimization of parameters 'a' and 'b' in the DCC-GJRGARCH model, resulting in significant performance inefficiencies. This update ensures that all necessary conditions are verified before proceeding with calculations, thus enhancing the model's efficiency and reliability. --- src/frds/algorithms/_garch.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/frds/algorithms/_garch.py b/src/frds/algorithms/_garch.py index 9857305..2f33111 100644 --- a/src/frds/algorithms/_garch.py +++ b/src/frds/algorithms/_garch.py @@ -359,6 +359,9 @@ def fit(self) -> Parameters: Returns: List[float]: [mu, omega, alpha, gamma, beta, loglikelihood] """ + # No repeated estimation? + if self.estimation_success: + return starting_vals = self.preparation() bounds = [ From af82225e49e575dc88270aab971221bd0d0570c6 Mon Sep 17 00:00:00 2001 From: lianghx123 <68636898+lianghx123@users.noreply.github.com> Date: Tue, 2 Jul 2024 20:27:11 +0800 Subject: [PATCH 3/5] Improved performance Vectorize the grid search method for finding initialization parameters and use the Numba JIT function to wrap the volatility model for improved performance. --- src/frds/algorithms/_garch.py | 220 ++++++++++++++++++++++++++-------- 1 file changed, 170 insertions(+), 50 deletions(-) diff --git a/src/frds/algorithms/_garch.py b/src/frds/algorithms/_garch.py index 2f33111..ab2dbb5 100644 --- a/src/frds/algorithms/_garch.py +++ b/src/frds/algorithms/_garch.py @@ -1,6 +1,6 @@ import warnings import itertools -from typing import List +from typing import List, Union from dataclasses import dataclass import numpy as np from scipy.optimize import minimize, OptimizeResult @@ -147,10 +147,10 @@ def loglikelihood_model( """ if zero_mean: resids = self.returns - self.sigma2 = self.compute_variance(params, resids, backcast, var_bounds) + self.sigma2 = compute_variance(params, resids, backcast, var_bounds) else: resids = self.returns - params[0] # params[0] is mu - self.sigma2 = self.compute_variance( + self.sigma2 = compute_variance( params[1:], resids, backcast, var_bounds ) return -self.loglikelihood(resids, self.sigma2) @@ -168,7 +168,10 @@ def loglikelihood(resids: np.ndarray, sigma2: np.ndarray) -> float: float: log-likelihood """ l = -0.5 * (np.log(2 * np.pi) + np.log(sigma2) + resids**2.0 / sigma2) - return np.sum(l) + if len(sigma2.shape)>1: + return np.sum(l, axis=1) + else: + return np.sum(l) @staticmethod def compute_variance( @@ -209,28 +212,58 @@ def starting_values(self, resids: np.ndarray) -> List[float]: Returns: List[float]: [omega, alpha, beta] """ - # Candidate target persistence - we don't have a prior knowledge here - persistence_grid = [0.5, 0.7, 0.9, 0.98] - # Candidate alpha values - alpha_grid = [0.01, 0.05, 0.1, 0.2] - # Sample variance - var = self.returns.var() - # Backcast initial value for the volaility process - initial_params = [] - max_likelihood = -np.inf - for alpha, p in itertools.product(*[alpha_grid, persistence_grid]): - # Note that (alpha+beta) is the "persistence" - beta = p - alpha - # Note that the unconditional variance is omega/(1-persistence). - # As we are guessing the initial value, we use the sample variance - # as the unconditional variance to guess the omega value. - omega = var * (1 - p) - params = [omega, alpha, beta] - sigma2 = self.compute_variance( - params, resids, self.backcast_value, self.var_bounds - ) - if -self.loglikelihood(resids, sigma2) > max_likelihood: - initial_params = params + def hierarchical_grid_search(resids, initial_range, base_grid_size=10, depth=2): + current_range = initial_range + best_params = None + + for i in range(depth): + if i != 0: + # 更新搜索范围为最佳参数周围的较小区域 + search_width = (current_range[:, 1] - current_range[:, 0]) / base_grid_size + current_range = np.array([ + [best_params[1] - search_width[0]*2, best_params[1] + search_width[0]*2], + [best_params[2] - search_width[1]*2, best_params[2] + search_width[1]*2] + ]) + # 创建当前深度的网格 + alpha_grid = np.linspace(current_range[0, 0], current_range[0, 1], base_grid_size) + p_grid = np.linspace(current_range[1, 0], current_range[1, 1], base_grid_size) + alpha_values, p_values = np.meshgrid(alpha_grid, p_grid) + parameters = np.stack([alpha_values.ravel(), p_values.ravel()], axis=1) + + # 过滤无效组合 + valid_indices = np.logical_and(parameters[:, 1] - parameters[:, 0] > 0, # beta > 0 + parameters[:, 0] > 0) # alpha > 0 + valid_parameters = parameters[valid_indices] + + # 得到三个参数组成的二维矩阵 + alpha = valid_parameters[:, 0] + p = valid_parameters[:, 1] + beta = p - alpha + omega = self.returns.var() * (1 - p) + + omega = omega.reshape(-1, 1) + alpha = alpha.reshape(-1, 1) + beta = beta.reshape(-1, 1) + + # Combine them into a new matrix + valid_parameters = np.hstack((omega, alpha, beta)) + sigma2 = compute_variance( + valid_parameters, resids, self.backcast_value, self.var_bounds + ) + ll_values = self.loglikelihood(resids, sigma2) + max_ll_index = np.argmax(ll_values) + best_params = valid_parameters[max_ll_index] + + return best_params + + # 初始化参数范围 + initial_range = np.array([[0.01, 0.2], [0.5, 0.99]]) # [a_range, b_range] + base_grid_size = 5 # 初始网格大小 + depth = 1 # 搜索深度，即细化的次数 + + # 调用函数 + best_parameters = hierarchical_grid_search(resids, initial_range, base_grid_size, depth) + initial_params = best_parameters return initial_params @staticmethod @@ -474,32 +507,119 @@ def starting_values(self, resids: np.ndarray) -> List[float]: Returns: List[float]: [omega, alpha, gamma, beta] """ - # Candidate target persistence - we don't have a prior knowledge here - persistence_grid = [0.5, 0.7, 0.9, 0.98] - # Candidate alpha values - alpha_grid = [0.01, 0.05, 0.1, 0.2] - gamma_grid = alpha_grid - # Sample variance - var = self.returns.var() - # Backcast initial value for the volaility process - var_bounds = self.variance_bounds(resids) - backcast = self.backcast(resids) - initial_params = [] - max_likelihood = -np.inf - # fmt: off - for alpha, gamma, p in itertools.product(*[alpha_grid, gamma_grid, persistence_grid]): - # Note that (alpha+beta) is the "persistence" - beta = p - alpha - gamma/2 - # Note that the unconditional variance is omega/(1-persistence). - # As we are guessing the initial value, we use the sample variance - # as the unconditional variance to guess the omega value. - omega = var * (1 - p) - params = [omega, alpha, gamma, beta] - sigma2 = self.compute_variance(params, resids, backcast, var_bounds) - if -self.loglikelihood(resids, sigma2) > max_likelihood: - initial_params = params + def hierarchical_grid_search(resids, initial_range, base_grid_size=10, depth=2): + current_range = initial_range + best_params = None + + for i in range(depth): + if i != 0: + # 更新搜索范围为最佳参数周围的较小区域 + search_width = (current_range[:, 1] - current_range[:, 0]) / base_grid_size + current_range = np.array([ + [best_params[1] - search_width[0]*2, best_params[1] + search_width[0]*2], + [best_params[2] - search_width[1]*2, best_params[2] + search_width[1]*2], + [best_params[3] - search_width[2]*2, best_params[3] + search_width[2]*2] + ]) + # 创建当前深度的网格 + alpha_grid = np.linspace(current_range[0, 0], current_range[0, 1], base_grid_size) + gamma_grid = np.linspace(current_range[1, 0], current_range[1, 1], base_grid_size) + p_grid = np.linspace(current_range[2, 0], current_range[2, 1], base_grid_size) + alpha_values, gamma_values, p_values = np.meshgrid(alpha_grid, gamma_grid, p_grid) + parameters = np.stack([alpha_values.ravel(), gamma_values.ravel(), p_values.ravel()], axis=1) + # 过滤无效组合 + valid_indices = np.logical_and(parameters[:, 2] - parameters[:, 0] - parameters[:, 1]/2 > 0, # beta > 0 + parameters[:, 1] > 0, # gamma>0 + parameters[:, 0] > 0) # alpha > 0 + valid_parameters = parameters[valid_indices] + + # 得到三个参数组成的二维矩阵 + alpha = valid_parameters[:, 0] + gamma = valid_parameters[:, 1] + p = valid_parameters[:, 2] + beta = p - alpha - gamma/2 + omega = self.returns.var() * (1 - p) + + omega = omega.reshape(-1, 1) + alpha = alpha.reshape(-1, 1) + gamma = gamma.reshape(-1, 1) + beta = beta.reshape(-1, 1) + + # Combine them into a new matrix + valid_parameters = np.hstack((omega, alpha, gamma, beta)) + sigma2 = compute_variance( + valid_parameters, resids, self.backcast_value, self.var_bounds + ) + ll_values = self.loglikelihood(resids, sigma2) + max_ll_index = np.argmax(ll_values) + best_params = valid_parameters[max_ll_index] + + return best_params + + # 初始化参数范围 + initial_range = np.array([[0.01, 0.2], [0.01, 0.2], [0.5, 0.99]]) # [a_range, b_range] + base_grid_size = 5 # 初始网格大小 + depth = 1 # 搜索深度，即细化的次数 + + # 调用函数 + best_parameters = hierarchical_grid_search(resids, initial_range, base_grid_size, depth) + initial_params = best_parameters return initial_params +from numba import jit +@jit(nopython=True, cache=True) +def bounds_check(sigma2: np.float64|np.ndarray, var_bounds:np.ndarray) -> np.float64|np.ndarray: + lower, upper = var_bounds[0], var_bounds[1] + if isinstance(sigma2, np.float64): + if sigma2 < lower: + sigma2 = lower + if sigma2 > upper: + if not np.isinf(sigma2): + sigma2 = upper + np.log(sigma2 / upper) + else: + sigma2 = upper + 1000.0 + else: + sigma2 = np.maximum(lower, sigma2) + sigma2 = np.where(sigma2 > upper, + np.where(np.isinf(sigma2), upper + 1000.0, upper + np.log(sigma2 / upper)), + sigma2) + return sigma2 + +@jit(nopython=True, cache=True) +def compute_variance(params: np.ndarray, resids: np.ndarray, backcast: np.float64, var_bounds: np.ndarray) -> np.ndarray: + if params.shape[-1] > 3: + gjrgarch = True + else: + gjrgarch = False + if params.ndim > 1: + omega = params[:, 0] + alpha = params[:, 1] + if gjrgarch: + gamma = params[:, 2] + beta = params[:, 3] + else: + gamma = np.zeros(omega.shape) + beta = params[:, 2] + sigma2 = np.zeros((omega.shape[0], resids.shape[0])) + sigma2[:, 0] = omega + (alpha + gamma/2 + beta) * backcast + resids_squared = resids ** 2 + for t in range(1, resids.shape[0]): + sigma2[:, t] = omega + alpha * resids_squared[t - 1] + beta * sigma2[:, t - 1] + sigma2[:, t] += gamma * resids_squared[t - 1] if resids[t - 1] < 0 else np.zeros_like(omega, dtype=np.float64) + sigma2[:, t] = bounds_check(sigma2[:, t], var_bounds[t]) + else: + if gjrgarch: + omega, alpha, gamma, beta = params + else: + omega, alpha, beta = params + gamma = 0 + sigma2 = np.zeros(resids.shape) + sigma2[0] = omega + (alpha + gamma/2 + beta) * backcast + for t in range(1, resids.shape[0]): + sigma2[t] = omega + alpha * (resids[t - 1] ** 2) + beta * sigma2[t - 1] + sigma2[t] += gamma * resids[t - 1] ** 2 if resids[t - 1] < 0 else 0 + sigma2[t] = bounds_check(sigma2[t], var_bounds[t]) + return sigma2 + if __name__ == "__main__": import pandas as pd From 2e815ecfdb156ae296e88c456c073e3b767de513 Mon Sep 17 00:00:00 2001 From: lianghx123 <68636898+lianghx123@users.noreply.github.com> Date: Tue, 2 Jul 2024 20:34:31 +0800 Subject: [PATCH 4/5] Improved performance Vectorize the grid search method for finding initialization parameters and use the Numba JIT function to wrap the dcc model for improved performance. --- src/frds/algorithms/_mgarch.py | 117 +++++++++++++++++++++++++++++---- 1 file changed, 103 insertions(+), 14 deletions(-) diff --git a/src/frds/algorithms/_mgarch.py b/src/frds/algorithms/_mgarch.py index f38fdf3..b5d40f4 100644 --- a/src/frds/algorithms/_mgarch.py +++ b/src/frds/algorithms/_mgarch.py @@ -304,17 +304,44 @@ def starting_values(self) -> Tuple[float, float]: Returns: Tuple[float, float]: [a, b] """ - a_grid = np.linspace(0.01, 0.5, 10) - b_grid = np.linspace(0.6, 0.95, 10) - max_ll = -np.inf - initial_values = [0.01, 0.8] - for a, b in itertools.product(a_grid, b_grid): - if a + b >= 1: - continue - ll = -self.loglikelihood_model([a, b]) - if ll > max_ll: - max_ll = ll - initial_values = [a, b] + def hierarchical_grid_search(loglikelihood_model, initial_range, base_grid_size=100, depth=2): + current_range = initial_range + best_params = None + + for i in range(depth): + # 创建当前深度的网格 + a_grid = np.linspace(current_range[0, 0], current_range[0, 1], base_grid_size) + b_grid = np.linspace(current_range[1, 0], current_range[1, 1], base_grid_size) + a_values, b_values = np.meshgrid(a_grid, b_grid) + parameters = np.stack([a_values.ravel(), b_values.ravel()], axis=1) + + # 过滤无效组合 + valid_indices = np.logical_and(parameters[:, 0] + parameters[:, 1] < 1, # a + b < 1 + parameters[:, 0] >= 0.001, # a >= 0.001 + parameters[:, 1] >= 0.001) # b >= 0.001 + valid_parameters = parameters[valid_indices] + + # 计算对数似然值 + ll_values = -loglikelihood_model(valid_parameters) + max_ll_index = np.argmax(ll_values) + best_params = valid_parameters[max_ll_index] + + # 更新搜索范围为最佳参数周围的较小区域 + search_width = (current_range[:, 1] - current_range[:, 0]) / base_grid_size + current_range = np.array([ + [best_params[0] - search_width[0]*2, best_params[0] + search_width[0]*2], + [best_params[1] - search_width[1]*2, best_params[1] + search_width[1]*2] + ]) + return best_params + + # 初始化参数范围 + initial_range = np.array([[0.001, 0.5], [0.6, 1]]) # [a_range, b_range] + base_grid_size = 10 # 初始网格大小 + depth = 1 # 搜索深度，即细化的次数 + + # 调用函数 + best_parameters = hierarchical_grid_search(self.loglikelihood_model, initial_range, base_grid_size, depth) + initial_values = best_parameters return initial_values def loglikelihood_model(self, params: np.ndarray) -> float: @@ -327,7 +354,12 @@ def loglikelihood_model(self, params: np.ndarray) -> float: float: negative log-likelihood """ # fmt: off - a, b = params + if len(params.shape)>1: + a = params[:, 0] # 第一列 + b = params[:, 1] # 第二列 + else: + a = params[0] + b = params[1] resids1 = self.model1.resids resids2 = self.model2.resids sigma2_1 = self.model1.sigma2 @@ -340,7 +372,13 @@ def loglikelihood_model(self, params: np.ndarray) -> float: l1 = self.model1.parameters.loglikelihood + self.model2.parameters.loglikelihood # The loglikelihood of the correlation component (Step 2) - rho = self.conditional_correlations(a, b) + self.model1.fit() # in case it was not estimated, no performance loss + self.model2.fit() # in case it was not estimated + resids1 = self.model1.resids + resids2 = self.model2.resids + sigma2_1 = self.model1.sigma2 + sigma2_2 = self.model2.sigma2 + rho = conditional_correlations(resids1, resids2, sigma2_1, sigma2_2, a, b) # TODO: rare case rho is out of bounds rho = np.clip(rho, -0.99, 0.99) @@ -349,7 +387,10 @@ def loglikelihood_model(self, params: np.ndarray) -> float: + np.log(1 - rho ** 2) + (z1 ** 2 + z2 ** 2 - 2 * rho * z1 * z2) / (1 - rho ** 2) ) - l2 = np.sum(log_likelihood_terms) + if len(params.shape)>1: + l2 = np.sum(log_likelihood_terms, axis=1) + else: + l2 = np.sum(log_likelihood_terms) negative_loglikelihood = - (l1 + l2) return negative_loglikelihood @@ -473,6 +514,54 @@ def fit(self) -> Parameters: """ return super().fit() +from numba import jit +@jit(nopython=True, cache=True) +def conditional_correlations(resids1, resids2, sigma2_1, sigma2_2, a, b) -> np.ndarray: + # z1 and z2 are standardized residuals + z1 = resids1 / np.sqrt(sigma2_1) + z2 = resids2 / np.sqrt(sigma2_2) + Q_bar = np.corrcoef(z1, z2) + q_11_bar, q_12_bar, q_22_bar = Q_bar[0, 0], Q_bar[0, 1], Q_bar[1, 1] + T = z1.shape[0] + z1_squared = np.square(z1) + z2_squared = np.square(z2) + z1_z2 = np.multiply(z1, z2) + coeff = 1 - a - b + + if isinstance(a, float): + q11 = np.empty_like(z1) + q12 = np.empty_like(z1) + q22 = np.empty_like(z1) + rho = np.zeros_like(z1) + q11[0] = q_11_bar + q22[0] = q_22_bar + q12[0] = q_12_bar + rho[0] = q12[0] / np.sqrt(q11[0] * q22[0]) + + for t in range(1, T): + q11[t] = coeff * q_11_bar + a * z1_squared[t - 1] + b * q11[t - 1] + q22[t] = coeff * q_22_bar + a * z2_squared[t - 1] + b * q22[t - 1] + q12[t] = coeff * q_12_bar + a * z1_z2[t - 1] + b * q12[t - 1] + rho[1:] = q12[1:] / np.sqrt(q11[1:] * q22[1:]) + + else: + q11 = np.empty((a.shape[0], z1.shape[0])) + q12 = np.empty((a.shape[0], z1.shape[0])) + q22 = np.empty((a.shape[0], z1.shape[0])) + rho = np.zeros((a.shape[0], z1.shape[0])) + q11[:, 0] = q_11_bar + q22[:, 0] = q_22_bar + q12[:, 0] = q_12_bar + rho[:, 0] = q12[:, 0] / np.sqrt(q11[:, 0] * q22[:, 0]) + + for t in range(1, T): + q11[:, t] = coeff * q_11_bar + a * z1_squared[t - 1] + b * q11[:, t - 1] + q22[:, t] = coeff * q_22_bar + a * z2_squared[t - 1] + b * q22[:, t - 1] + q12[:, t] = coeff * q_12_bar + a * z1_z2[t - 1] + b * q12[:, t - 1] + rho[:, 1:] = q12[:, 1:] / np.sqrt(q11[:, 1:] * q22[:, 1:]) + + return rho + if __name__ == "__main__": import pandas as pd From d6ff5871e68b10a3f32fb9711189ebde59e0261f Mon Sep 17 00:00:00 2001 From: lianghx123 <68636898+lianghx123@users.noreply.github.com> Date: Tue, 2 Jul 2024 20:40:31 +0800 Subject: [PATCH 5/5] Enhanced performance Vectorize the S-sample calculations and computation process of LRMES to enhance performance. --- src/frds/measures/_long_run_mes.py | 155 ++++++++++++++--------------- 1 file changed, 72 insertions(+), 83 deletions(-) diff --git a/src/frds/measures/_long_run_mes.py b/src/frds/measures/_long_run_mes.py index 3690679..5a8eedc 100644 --- a/src/frds/measures/_long_run_mes.py +++ b/src/frds/measures/_long_run_mes.py @@ -1,6 +1,7 @@ from typing import Tuple import numpy as np from frds.algorithms import GJRGARCHModel, GJRGARCHModel_DCC +from frds.algorithms._mgarch import conditional_correlations USE_CPP_EXTENSION = True try: @@ -62,7 +63,13 @@ def estimate(self, h=22, S=10_000, C=-0.1, random_seed=42) -> float: market_resids = self.market_model.resids # Conditional correlations a, b = self.dcc_model.parameters.a, self.dcc_model.parameters.b - rho = self.dcc_model.conditional_correlations(a, b) + self.dcc_model.model1.fit() # in case it was not estimated, no performance loss + self.dcc_model.model2.fit() # in case it was not estimated + resids1 = self.dcc_model.model1.resids + resids2 = self.dcc_model.model2.resids + sigma2_1 = self.dcc_model.model1.sigma2 + sigma2_2 = self.dcc_model.model2.sigma2 + rho = conditional_correlations(resids1, resids2, sigma2_1, sigma2_2, a, b) # Standarized residuals z_m = market_resids / np.sqrt(market_variances) z_i = firm_resids / np.sqrt(firm_variances) @@ -73,65 +80,32 @@ def estimate(self, h=22, S=10_000, C=-0.1, random_seed=42) -> float: # Sample with replacement S*h innovations # sample = sample.T[rng.choice(sample.shape[1], (S, h), replace=True)] + idx = rng.choice(len(xi_i), size=(S, h), replace=True) + # 使用高级索引和广播来生成 sample 数组 sample = np.zeros((S, h, 2)) - # Sample with replacement - for s in range(S): - idx = rng.choice(len(xi_i), size=h, replace=True) - sample[s, :, 0] = xi_i[idx] - sample[s, :, 1] = z_m[idx] + sample[:, :, 0] = xi_i[idx] + sample[:, :, 1] = z_m[idx] assert sample.shape == (S, h, 2) Q_bar = np.corrcoef(z_i, z_m) - firm_avg_return = 0.0 - n_systemic_event = 0 - for s in range(S): - # Each simulation - inv = sample[s, :, :] # shape=(h,2) - assert inv.shape == (h, 2) - if USE_CPP_EXTENSION: - firm_return, systemic_event = ext.simulation( - inv, - C, - self.firm_model.sigma2[-1], - self.market_model.sigma2[-1], - self.firm_model.resids[-1], - self.market_model.resids[-1], - self.dcc_model.parameters.a, - self.dcc_model.parameters.b, - rho[-1], - Q_bar, - self.firm_model.parameters.mu, - self.firm_model.parameters.omega, - self.firm_model.parameters.alpha, - self.firm_model.parameters.gamma, - self.firm_model.parameters.beta, - self.market_model.parameters.mu, - self.market_model.parameters.omega, - self.market_model.parameters.alpha, - self.market_model.parameters.gamma, - self.market_model.parameters.beta, - ) - else: - firm_return, systemic_event = self.simulation( - innovation=inv, - C=C, - firm_var=self.firm_model.sigma2[-1], - mkt_var=self.market_model.sigma2[-1], - firm_resid=self.firm_model.resids[-1], - mkt_resid=self.market_model.resids[-1], - a=self.dcc_model.parameters.a, - b=self.dcc_model.parameters.b, - rho=rho[-1], - Q_bar=Q_bar, - ) - firm_avg_return += firm_return - n_systemic_event += systemic_event - - if n_systemic_event == 0: - return 0.0 - return -firm_avg_return / n_systemic_event + systemic_event_firm_return = self.simulation( + innovation=sample, + C=C, + firm_var=self.firm_model.sigma2[-1], + mkt_var=self.market_model.sigma2[-1], + firm_resid=self.firm_model.resids[-1], + mkt_resid=self.market_model.resids[-1], + a=self.dcc_model.parameters.a, + b=self.dcc_model.parameters.b, + rho=rho[-1], + Q_bar=Q_bar, + ) + systemic_event_firm_return = systemic_event_firm_return[systemic_event_firm_return != False] + + # 计算这些值的均值 + return -np.mean(systemic_event_firm_return) def simulation( self, @@ -167,48 +141,69 @@ def simulation( q_i_bar = Q_bar[0, 0] q_m_bar = Q_bar[1, 1] q_im_bar = Q_bar[1, 0] - q_i, q_m, q_im = 1.0, 1.0, rho + S, h, d = innovation.shape + q_i = np.full(S, 1) + q_m = np.full(S, 1) + q_im = np.full(S, rho) pi = self.firm_model.parameters mu_i = pi.mu - omega_i, alpha_i, gamma_i, beta_i = pi.omega, pi.alpha, pi.gamma, pi.beta + + # omega_i, alpha_i, gamma_i, beta_i = pi.omega, pi.alpha, pi.gamma, pi.beta + omega_i = np.full(S, pi.omega) + alpha_i = np.full(S, pi.alpha) + gamma_i = np.full(S, pi.gamma) + beta_i = np.full(S, pi.beta) pm = self.market_model.parameters mu_m = pm.mu - omega_m, alpha_m, gamma_m, beta_m = pm.omega, pm.alpha, pm.gamma, pm.beta + # omega_m, alpha_m, gamma_m, beta_m = pm.omega, pm.alpha, pm.gamma, pm.beta + omega_m = np.full(S, pm.omega) + alpha_m = np.full(S, pm.alpha) + gamma_m = np.full(S, pm.gamma) + beta_m = np.full(S, pm.beta) - firm_return = np.empty(shape=len(innovation)) - mkt_return = np.empty(shape=len(innovation)) + firm_return = np.empty((S, h)) + mkt_return = np.empty((S, h)) # lagged residuals - resid_i_tm1 = firm_resid - resid_m_tm1 = mkt_resid + # resid_i_tm1 = firm_resid + # resid_m_tm1 = mkt_resid + resid_i_tm1 = np.full(S, firm_resid) + resid_m_tm1 = np.full(S, mkt_resid) # lagged conditonal variance - firm_var_tm1 = firm_var - mkt_var_tm1 = mkt_var + # firm_var_tm1 = firm_var + # mkt_var_tm1 = mkt_var + firm_var_tm1 = np.full(S, firm_var) + mkt_var_tm1 = np.full(S, mkt_var) # lagged standarized residuals - firm_innov_tm1 = resid_i_tm1 / np.sqrt(firm_var_tm1) - mkt_innov_tm1 = resid_m_tm1 / np.sqrt(mkt_var_tm1) + # firm_innov_tm1 = resid_i_tm1 / np.sqrt(firm_var_tm1) + # mkt_innov_tm1 = resid_m_tm1 / np.sqrt(mkt_var_tm1) + firm_innov_tm1 = np.full(S, firm_resid / np.sqrt(firm_var)) + mkt_innov_tm1 = np.full(S, mkt_resid / np.sqrt(mkt_var)) + coff = (1 - a - b) - for h in range(len(innovation)): + for h in range(innovation.shape[-2]): # fmt: off # Each iteration is a one-step-ahead forecast # conditional variance this time firm_var_t = omega_i + alpha_i * (resid_i_tm1**2) + beta_i * firm_var_tm1 - firm_var_t += gamma_i * (resid_i_tm1**2) if resid_i_tm1 < 0 else 0 + # firm_var_t += gamma_i * (resid_i_tm1**2) if resid_i_tm1 < 0 else 0 + firm_var_t += np.where(resid_i_tm1 < 0, gamma_i * (resid_i_tm1**2), 0) mkt_var_t = omega_m + alpha_m * (resid_m_tm1**2) + beta_m * mkt_var_tm1 - mkt_var_t += gamma_m * (resid_m_tm1**2) if resid_m_tm1 < 0 else 0 + # mkt_var_t += gamma_m * (resid_m_tm1**2) if resid_m_tm1 < 0 else 0 + mkt_var_t += np.where(resid_m_tm1 < 0, gamma_m * (resid_m_tm1**2), 0) # conditional correlation this time - q_i = (1 - a - b) * q_i_bar + a * firm_innov_tm1**2 + b * q_i - q_m = (1 - a - b) * q_m_bar + a * mkt_innov_tm1**2 + b * q_m - q_im = (1 - a - b) * q_im_bar + a * firm_innov_tm1*mkt_innov_tm1 + b * q_im + q_i = coff * q_i_bar + a * firm_innov_tm1**2 + b * q_i + q_m = coff * q_m_bar + a * mkt_innov_tm1**2 + b * q_m + q_im = coff * q_im_bar + a * firm_innov_tm1*mkt_innov_tm1 + b * q_im rho_h = q_im / np.sqrt(q_i * q_m) # innovations this time - firm_innov = innovation[h, 0] - mkt_innov = innovation[h, 1] + firm_innov = innovation[:, h, 0] + mkt_innov = innovation[:, h, 1] # market excess return # or, residual this time, conditional volatility * innovation (z_m) @@ -221,8 +216,8 @@ def simulation( rho_h * mkt_innov + np.sqrt(1 - rho_h**2) * firm_innov ) - mkt_return[h] = mu_m + epsilon_m - firm_return[h] = mu_i + epsilon_i + mkt_return[:, h] = mu_m + epsilon_m + firm_return[:, h] = mu_i + epsilon_i firm_var_tm1 = firm_var_t mkt_var_tm1 = mkt_var_t @@ -237,12 +232,6 @@ def simulation( # systemic event if over the prediction horizon, # the market falls by more than C - systemic_event = np.exp(np.sum(mkt_return)) - 1 < C + systemic_event = np.exp(np.sum(mkt_return, axis=1)) - 1 < C # no need to simulate firm returns if there is no systemic event - if not systemic_event: - return 0.0, False - - # firm return over the horizon - firmret = np.exp(np.sum(firm_return)) - 1 - - return firmret, True + return np.where(systemic_event, np.exp(np.sum(firm_return, axis=1)) - 1, False)