WIP: EVT VaR and ES.

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()

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,22 +121,19 @@ 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.
         """
         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

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.

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):

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",
                                 )

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()

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()