Working Vine Copula modeling.

Sinbad-The-Sailor · Aug 6, 2023 · de47bb5 · de47bb5
1 parent 07291ba
commit de47bb5
Show file tree

Hide file tree

Showing 8 changed files with 141 additions and 59 deletions.
diff --git a/.todo b/.todo
@@ -1,6 +1,6 @@
 # 1. For AR(p) models the number of observations, n, is assumed to be sufficiently
 #    large to use the conditional log-likelihood instead of full likelihood.
 
-# 2. For more advanced processes, especially dependent on graident based maximum
-#    likelihood an integration with pytorch has to be made to utilize automatic
-#    differentiation.
+# 2. For GARCH models consider estimating the mean, mu, too.
+
+# 3. Clean up and divide torch tensors and numpy arrays appropriately.
diff --git a/src2/autograd.ipynb → src2/examples/autograd.ipynb b/src2/autograd.ipynb → src2/examples/autograd.ipynb
diff --git a/src2/data.csv → src2/examples/data.csv b/src2/data.csv → src2/examples/data.csv
diff --git a/src2/run.ipynb → src2/examples/run.ipynb b/src2/run.ipynb → src2/examples/run.ipynb
diff --git a/src2/testing.ipynb → src2/examples/testing.ipynb b/src2/testing.ipynb → src2/examples/testing.ipynb
diff --git a/src2/main.py b/src2/main.py
@@ -1,68 +1,72 @@
 # -*- coding: utf-8 -*-
 import torch
-import pandas as pd
 import numpy as np
 
 from datetime import datetime
 from matplotlib import pyplot as plt
 
 from utils.stock_factory import StockFactory
-from utils.instruments import Stock
-from utils.config import STOCK_ADMISSIBLE_MODELS
+from simulator.simulator import Simulator
 
-from models.ar import AR
-from models.garch import GARCH
 
 
-
-start = datetime.strptime("2005-05-01", r"%Y-%m-%d")
+start = datetime.strptime("2013-05-01", r"%Y-%m-%d")
 end = datetime.strptime("2023-06-01", r"%Y-%m-%d")
-instrument_specification = ["XOM", "GS", "T"]
+instrument_specification = ("XOM",
+                            "GS",
+                            "T",
+                            "AAPL",
+                            "MSFT",
+                            "PG",
+                            "K",
+                            "ADI",
+                            "GE",
+                            "AIG",
+                            "KO",
+                            "NKE",
+                            "BAC",
+                            "MMM",
+                            "AXP",
+                            "AMZN",
+                            )
 instrument_factory = StockFactory(tickers=instrument_specification,
                                   start=start,
                                   end=end)
 stocks = instrument_factory.build_stocks()
 
 
+number_of_simulations = 10
+time_steps = 100
 
-class Simulator:
-
-    def __init__(self, instruments: list):
-        self._instruments = instruments
-        self._calibrated = False
-
-    def calibrate(self):
-        self._calibrate_instruments()
-        self._calibrate_copula()
-        self._calibrated = True
-
-    def _calibrate_instruments(self):
-        for instrument in self._instruments:
-            if isinstance(instrument, Stock):
-                self._calibrate_stock(instrument)
+sim = Simulator(stocks)
+sim.calibrate()
+simulation_tensor = sim.run_simulation(time_steps, number_of_simulations)
 
-    def _calibrate_copula(self):
-        ...
 
-    def _calibrate_stock(self, stock):
-        current_aic = np.inf
-        risk_factor = stock.risk_factors[0]
-        data = risk_factor.price_history.log_returns
 
-        for model_name in STOCK_ADMISSIBLE_MODELS:
-            if model_name == "AR":
-                model = AR(data)
-                model.calibrate()
-                if model.aic < current_aic:
-                    risk_factor.model = model
+fig, ax = plt.subplots(4, 4)
+ix, iy = 0, 0
 
-            elif model_name == "GARCH":
-                model = GARCH(data)
-                model.calibrate()
-                if model.aic < current_aic:
-                    risk_factor.model = model
+for i, stock in enumerate(stocks):
+    stock_name = stock.identifier
+    past_prices = stock.risk_factors[0].price_history.mid_history[-500:]
+    past_time = range(len(past_prices))
+    ax[ix, iy].plot(past_time, past_prices)
+    ax[ix, iy].set_title(stock_name)
 
+    future_time = range(len(past_time) - 1, len(past_time) + time_steps)
+    inital_price = torch.tensor(past_prices[-1])
+    for k in range(number_of_simulations):
+        all_prices = torch.empty(time_steps+1)
+        prices = torch.exp(torch.cumsum(simulation_tensor[i,:,k], dim=0)) * inital_price
+        all_prices[0] = inital_price
+        all_prices[1:] = prices
+        ax[ix, iy].plot(future_time, all_prices, color="grey")
 
-    def run_simulation(time_steps: int) -> torch.Tensor:
-        # Check for succesful calibration, Throw an error otherwise.
-        ...
+    if ix != 3:
+        ix += 1
+    else:
+        iy += 1
+        ix = 0
+plt.tight_layout()
+plt.show()
diff --git a/src2/simulator/simulator.py b/src2/simulator/simulator.py
@@ -1,25 +1,87 @@
 # -*- coding: utf-8 -*-
 
 import torch
-import pandas as pd
+import numpy as np
+import pyvinecopulib as pv
 
-from utils.instruments import Asset
+from utils.config import VINE_COPULA_FAMILIES, VINE_COPULA_NUMBER_OF_THREADS, STOCK_ADMISSIBLE_MODELS
+from utils.instruments import Stock, RiskFactor
+from models.ar import AR
+from models.garch import GARCH
 
 class Simulator:
 
-    def __init__(self, instruments: list[Asset]):
+    def __init__(self, instruments: list):
         self._instruments = instruments
-        self.is_calibrated = False
+        self._calibrated = False
 
+    @property
+    def _uniform_samples(self) -> np.array:
+        samples = []
+        for instrument in self._instruments:
+            risk_factor = instrument.risk_factors[0]
+            model = risk_factor.model
+            samples.append(model.transform_to_uniform())
+        return np.stack(samples).T
 
+    @property
+    def _number_of_risk_factors(self) -> int:
+        return len(self._risk_factors)
+
+    @property
+    def _risk_factors(self) -> list[RiskFactor]:
+        risk_factors = []
+        for instrument in self._instruments:
+            risk_factors += instrument.risk_factors
+        return risk_factors
 
     def calibrate(self):
-        # Find models for each risk factor.
-        # Use a factory design pattern.
-        # If none are succesful throw an error.
-        ...
+        self._calibrate_instruments()
+        self._calibrate_copula()
+        self._calibrated = True
+
+    def run_simulation(self, time_steps: int, number_of_simulations: int) -> torch.Tensor:
+        number_of_risk_factors = self._number_of_risk_factors
+        size = (number_of_risk_factors, time_steps, number_of_simulations)
+        simulation_tensor = torch.empty(size=size)
+
+        print("number of simulations", number_of_simulations)
+        print("time steps", time_steps)
+        print("number of risk factors", self._number_of_risk_factors)
+        print("number of simulations from cop", len(self._coupla.simulate(time_steps)))
+
+        for n in range(number_of_simulations):
+            simulations = self._coupla.simulate(time_steps).T
+            for i, simulation in enumerate(simulations):
+                simulation_tensor[i,:,n] = self._risk_factors[i].model.transform_to_true(torch.tensor(simulation))
+        return simulation_tensor
+
+    def _calibrate_instruments(self):
+        for instrument in self._instruments:
+            if isinstance(instrument, Stock):
+                self._calibrate_stock(instrument)
+
+    def _calibrate_copula(self):
+        uniforms = self._uniform_samples
+        controls = pv.FitControlsVinecop(family_set=VINE_COPULA_FAMILIES,
+                                         num_threads=VINE_COPULA_NUMBER_OF_THREADS)
+        copula = pv.Vinecop(uniforms, controls=controls)
+        self._coupla = copula
+
+    def _calibrate_stock(self, stock):
+        current_aic = np.inf
+        risk_factor = stock.risk_factors[0]
+        data = risk_factor.price_history.log_returns
 
+        for model_name in STOCK_ADMISSIBLE_MODELS:
+            if model_name == "AR":
+                model = AR(data)
+                model.calibrate()
+                if model.aic < current_aic:
+                    risk_factor.model = model
 
-    def run_simulation(time_steps: int) -> torch.Tensor:
-        # Check for succesful calibration, Throw an error otherwise.
-        ...
+            elif model_name == "GARCH":
+                model = GARCH(data)
+                model.calibrate()
+                if model.aic < current_aic:
+                    risk_factor.model = model
diff --git a/src2/utils/config.py b/src2/utils/config.py
@@ -1,4 +1,5 @@
 # -*- coding: utf-8 -*-
+import pyvinecopulib as pv
 
 # Years and Months
 MINIMUM_PRICE_HISTORY = (1, 0)
@@ -7,6 +8,21 @@
 MAXIMUM_AR_ORDER = 5
 INITIAL_VARIANCE_GARCH_OBSERVATIONS = 20
 INITIAL_GARCH_PARAMETERS = (0.05, 0.90)
+STOCK_ADMISSIBLE_MODELS = (
+    "AR",
+    "GARCH"
+)
+VINE_COPULA_FAMILIES = (
+    pv.BicopFamily.gaussian,
+    pv.BicopFamily.clayton,
+    pv.BicopFamily.frank,
+    pv.BicopFamily.gumbel,
+    pv.BicopFamily.student,
+    pv.BicopFamily.bb1,
+    pv.BicopFamily.bb6,
+    pv.BicopFamily.bb7,
+    pv.BicopFamily.bb8,
+)
 
-STOCK_ADMISSIBLE_MODELS = ("AR",
-                           "GARCH")
+# Computations
+VINE_COPULA_NUMBER_OF_THREADS = 6