From bd010d6a079e03b38bb8d960dfcf2c03d605a0aa Mon Sep 17 00:00:00 2001 From: Axel Nilsson Date: Sat, 19 Aug 2023 13:07:14 +0200 Subject: [PATCH] WIP: EVT VaR and ES. --- src2/main.py | 10 ++- src2/models/ar.py | 8 -- src2/models/model.py | 1 - src2/simulator/simulator.py | 1 - src2/utils/config.py | 28 ++++--- src2/utils/generalized_pareto.py | 26 ++++++ src2/utils/risk_assessor.py | 140 +++++++++++++++++++++++++++++-- 7 files changed, 182 insertions(+), 32 deletions(-) diff --git a/src2/main.py b/src2/main.py index e89cd2d..e6659d3 100644 --- a/src2/main.py +++ b/src2/main.py @@ -33,7 +33,7 @@ portfolio = Portfolio(holdings, cash) # Simulation specifiation. -number_of_simulations = 1000 +number_of_simulations = 10000 time_steps = 5 # Simulator for returns and prices. @@ -46,9 +46,11 @@ # Risk Assesor of simulated portfolio. risk_assessor = RiskAssessor(portfolio, return_tensor, 2) -print(risk_assessor._weights) -print(risk_assessor._return_matrix) -print(risk_assessor._reduced_return_matrix) + +print(risk_assessor.value_at_risk(0.99)) +print(risk_assessor.expected_shortfall(0.99)) +print(risk_assessor.extreme_value_at_risk(0.99)) +print(risk_assessor.extreme_expected_shortfall(0.99)) import numpy as np plt.hist(np.array(risk_assessor._portfolio_return_scenarios)) plt.show() diff --git a/src2/models/ar.py b/src2/models/ar.py index e105062..64fbecb 100644 --- a/src2/models/ar.py +++ b/src2/models/ar.py @@ -40,7 +40,6 @@ class AR(Model): def __init__(self, time_series: pd.Series): super().__init__(time_series) - @property def parameters(self) -> list[float]: self._check_calibration() @@ -73,7 +72,6 @@ def _log_likelihood(self) -> torch.Tensor: return - (self._number_of_observations - self._order) / 2 * (torch.log(pi) + torch.log(variance)) - 1 / (2 * variance) * torch.sum(squared_difference) - def calibrate(self): self._build_lag_order() self._build_data_matricies() @@ -83,7 +81,6 @@ def calibrate(self): self._calibrated = True - def transform_to_true(self, uniform_sample: torch.Tensor) -> torch.Tensor: self._check_calibration() @@ -103,14 +100,12 @@ def transform_to_true(self, uniform_sample: torch.Tensor) -> torch.Tensor: return simulated_values[self._order:] - def transform_to_uniform(self) -> torch.Tensor: self._check_calibration() residuals = self._build_residuals() return Normal(0, 1).cdf(residuals) - def _check_unit_roots(self): """ Checks for unit roots outside the unit circle. @@ -126,7 +121,6 @@ def _check_unit_roots(self): if np.abs(root) >= 1: raise StationarityError(f"Root {root} outside of complex unit circle.") - def _estimate_standard_deviation(self): """ Estimates unconditional standard deviation taking into account the order of the process. @@ -134,14 +128,12 @@ def _estimate_standard_deviation(self): variance = 1 / (self._number_of_observations - self._order) * torch.sum(torch.square(self._data - torch.mean(self._data))) self._sigma = torch.sqrt(variance) - def _solve_least_squares(self): solution = torch.linalg.solve(self._R, self._Q.T @ self._b) self._solution = solution self._mu = solution[0] self._phi = solution[1:] - def _build_data_matricies(self): """ The data matrix is build column-wise to represent regression matrix for least diff --git a/src2/models/model.py b/src2/models/model.py index 1a3795e..a91ed05 100644 --- a/src2/models/model.py +++ b/src2/models/model.py @@ -83,7 +83,6 @@ def bic(self) -> torch.Tensor: return number_of_parameters * torch.log(number_of_observations) - 2 * log_likelihood - def _check_calibration(self): """ Checks if successful calibration has been made. diff --git a/src2/simulator/simulator.py b/src2/simulator/simulator.py index 149bd07..fa05355 100644 --- a/src2/simulator/simulator.py +++ b/src2/simulator/simulator.py @@ -73,7 +73,6 @@ def run_simulation(self, time_steps: int, number_of_simulations: int) -> torch.T simulation_tensor[i,:,n] = self._risk_factors[i].model.transform_to_true(torch.tensor(simulation)) self._return_tensor = simulation_tensor - def _calibrate_instruments(self): for instrument in self._instruments: if isinstance(instrument, Stock): diff --git a/src2/utils/config.py b/src2/utils/config.py index 62f99a2..73400bc 100644 --- a/src2/utils/config.py +++ b/src2/utils/config.py @@ -27,6 +27,10 @@ # Computations VINE_COPULA_NUMBER_OF_THREADS = 6 +# Extreme Value Theory POT Threshold +EVT_THRESHOLD = 0.95 +GEV_INITIAL_SOLUTION = (0.1, 0.01) + # Optimizer DEFAULT_SOLVER = "ipopt" @@ -35,16 +39,16 @@ "GS", "T", "AAPL", - "MSFT", - "PG", - "K", - "ADI", - "GE", - "AIG", - "KO", - "NKE", - "BAC", - "MMM", - "AXP", - "AMZN", + # "MSFT", + # "PG", + # "K", + # "ADI", + # "GE", + # "AIG", + # "KO", + # "NKE", + # "BAC", + # "MMM", + # "AXP", + # "AMZN", ) diff --git a/src2/utils/generalized_pareto.py b/src2/utils/generalized_pareto.py index e69de29..6bee071 100644 --- a/src2/utils/generalized_pareto.py +++ b/src2/utils/generalized_pareto.py @@ -0,0 +1,26 @@ +# -*- coding: utf-8 -*- +import torch + + +class GeneralizedPareto: + """ + Representing all functionalities tied to a generalized pareto distribution. The object is instantiated in the + file, used by importing the file and using the instance of 'gp' ('generalized pareto'). + """ + + @staticmethod + def pdf(x: torch.Tensor, xi: torch.Tensor, beta: torch.Tensor) -> torch.Tensor: + if xi == 0: + return 1 / beta * torch.exp(-x / beta) + else: + return 1 / beta * (1 + xi / beta * x) ** (-1 / xi - 1) + + @staticmethod + def cdf(x: torch.Tensor, xi: torch.Tensor, beta: torch.Tensor) -> torch.Tensor: + if xi == 0: + return 1 - torch.exp(-x / beta) + else: + return 1 - (1 + xi / beta * x) ** (-1 / xi) + + +gp = GeneralizedPareto() diff --git a/src2/utils/risk_assessor.py b/src2/utils/risk_assessor.py index 13d9a91..8c89eb5 100644 --- a/src2/utils/risk_assessor.py +++ b/src2/utils/risk_assessor.py @@ -1,19 +1,43 @@ # -*- coding: utf-8 -*- import torch +import numpy as np +from scipy.optimize import minimize + +from utils.config import EVT_THRESHOLD, GEV_INITIAL_SOLUTION from utils.portfolio import Portfolio class RiskAssessor: + """ + Used for the risk assessment of portfolio returns given an empirical return distribution. Risk assessment is found + through standard measurements of Value at Risk and Expected Shortfall with user-specified confidence levels. Such + measurements are calculated with Extreme Value Theory (EVT). + + Portfolio returns are translated into portfolio losses in the class. + + The class is responsible for slicing the return tensor into the shape of the portfolio at the given time step. + """ def __init__(self, portfolio: Portfolio, return_tensor: torch.Tensor, time_step: int): self._portfolio = portfolio self._return_tensor = return_tensor self._time_step = time_step + self._calibrated = False self._check_time_step() + @property + def _evt_threshold(self) -> torch.Tensor: + return torch.quantile(self._portfolio_loss_scenarios, torch.tensor(EVT_THRESHOLD)) + + @property + def _excess_portfolio_losses(self) -> torch.Tensor: + large_losses = self._portfolio_loss_scenarios[self._portfolio_loss_scenarios > self._evt_threshold] + excess_losses = large_losses - self._evt_threshold + return excess_losses + @property def _return_matrix(self) -> torch.Tensor: return self._return_tensor[:,self._time_step,:] @@ -45,22 +69,126 @@ def _portfolio_return_scenarios(self) -> torch.Tensor: def _portfolio_loss_scenarios(self) -> torch.Tensor: return -self._portfolio_return_scenarios + @property + def _constraints(self) -> list[dict]: + constraints = [] + observations = self._excess_portfolio_losses + + for observation in observations: + constraints.append({"type": "ineq", "fun": lambda x: 1 + x[0] / x[1] * observation}) + + return constraints + def _check_time_step(self): if self._time_step > self._maxmimum_time_step: raise ValueError(f"Time step {self._time_step} is above maximum: {self._maxmimum_time_step}.") + def _calibrate(self): + solution = minimize(self._cost_function, + np.array(GEV_INITIAL_SOLUTION), + jac=True, + constraints=self._constraints) + self._solution = solution + self._parameters = solution.x + self._calibrated = True + + def _cost_function(self, parameters: np.array) -> tuple[float, 2]: + """ + Internal log-loss function. Note, this returns the log-loss with associated gradient. + + Args: + params: parameters for the generalized pareto distribution + data: observations above threshold determined in EVT. + + Returns: the log-loss value and the gradient of the generalized pareto distribution. + """ + parameters = torch.tensor(parameters, requires_grad=True) + log_loss = 0 + observations = self._excess_portfolio_losses + number_of_observations = len(observations) + + for observation in observations: + log_loss = log_loss + torch.log(1 + parameters[0] / parameters[1] * observation) + + log_loss = number_of_observations * torch.log(parameters[1]) + (1 + 1 / parameters[0]) * log_loss + + log_loss.backward() + return log_loss.data.cpu().numpy(), parameters.grad.data.cpu().numpy() + + def extreme_value_at_risk(self, confidence_level: float) -> float: - ... + """ + Calculates the Value at Risk of portfolio returns (either in percentage terms or in absolute value) as specified + in the instantiation of the RiskAssessor, based on EVT. + + Args: + quantile: the quantile of the Value at Risk measurement. Note, portfolio losses are specified! + + Returns: EVT based Value at Risk for a given confidence level. + """ + if not self._calibrated: + self._calibrate() + confidence_level = torch.tensor(confidence_level) + total_number_of_observations = self._portfolio_loss_scenarios.size(dim=0) + number_of_excess_observations = self._excess_portfolio_losses.size(dim=0) + threshold = self._evt_threshold + parameters = torch.tensor(self._parameters) + xi, beta = parameters[0], parameters[1] + print(xi, beta) + + + extreme_var = threshold + beta / xi * ((total_number_of_observations / number_of_excess_observations) * (1 - confidence_level) ** (-xi) - 1) + + return extreme_var.item() + def extreme_expected_shortfall(self, confidence_level: float) -> float: - ... + """ + Calculates the Expected Shortfall of portfolio returns (either in percentage terms or in absolute value) as + specified in the instantiation of the RiskAssessor, based on EVT. + + Args: + quantile: the quantile of the Value at Risk measurement. Note, portfolio losses are specified! + + Returns: EVT based Expected Shortfall for a given confidence level. + """ + if not self._calibrated: + self._calibrate() + extreme_var = torch.tensor(self.extreme_value_at_risk(confidence_level)) + threshold = self._evt_threshold + parameters = torch.tensor(self._parameters) + xi, beta = parameters[0], parameters[1] + + print(xi, beta) + + extreme_es = extreme_var / (1 - xi) + (beta - threshold * xi) /(1 - xi) + + return extreme_es.item() + def value_at_risk(self, confidence_level: float) -> float: + """ + Calculates the Value at Risk of portfolio returns (either in percentage terms or in absolute value) as specified + in the instantiation of the RiskAssessor. + + Args: + quantile: the quantile of the Value at Risk measurement. Note, portfolio losses are specified! + + Returns: Value at Risk for a given confidence level. + """ confidence_level = torch.tensor(confidence_level) - return torch.quantile(self._portfolio_loss_scenarios, confidence_level) + return torch.quantile(self._portfolio_loss_scenarios, confidence_level).item() def expected_shortfall(self, confidence_level: float) -> float: - confidence_level = torch.tensor(confidence_level) + """ + Calculates the Expected Shortfall of portfolio returns (either in percentage terms or in absolute value) as + specified in the instantiation of the RiskAssessor. + + Args: + quantile: the quantile of the Value at Risk measurement. Note, portfolio losses are specified! + + Returns: Expected Shortfall for a given confidence level. + """ var = self.value_at_risk(confidence_level) - losses_greater_than_var = torch.le(self._portfolio_loss_scenarios, var) - return torch.mean(losses_greater_than_var) + losses_greater_than_var = self._portfolio_loss_scenarios[self._portfolio_loss_scenarios > var] + return torch.mean(losses_greater_than_var).item()