Factored out discount() function.

Changed hedge calculation to use average of all available simulated
prices, which allows model time step to more flexible (step through
monthly periodisation in steps of e.g. 7 days).

Also, tried to make Wait multi-inherit from Settlement and Fixing,
but it doesn't work for a reason that eluded me (simulation requirements
somehow get messed up). So there are some comments about that.

Also changed plots to use percentiles, because it looks more intuitive.

quantdsl/__init__.py

@@ -1,6 +1,9 @@
 
 __version__ = '1.2.0dev0'
 
+# Todo: Option to switch between single- and double-sided deltas (currently always double-sided).
+# Todo: Tidy up how args are passed into evaluate(), it seems correct, but also a bit ad hoc.
+# Todo: Support names in expressions being resolved by evaluation args (e.g. like 'observation_date' but more general).
 # Todo: StockMarket element.
 # Todo: Better report object.
 # Todo: Warning when a built-in is being overridden by a user defined function. Or perhaps that would be useful? In any case, make sure this is handled properly (currently the user defined function will just be ignore?).

quantdsl/interfaces/calcandplot.py

@@ -15,6 +15,7 @@
 import numpy
 import six
 from eventsourcing.domain.model.events import subscribe, unsubscribe
+from numpy.lib.nanfunctions import nanpercentile
 
 from quantdsl.application.base import DEFAULT_MAX_DEPENDENCY_GRAPH_SIZE, QuantDslApplication
 from quantdsl.application.with_multithreading_and_python_objects import \
@@ -25,10 +26,11 @@
 from quantdsl.domain.model.call_result import CallResult, ResultValueComputed, make_call_result_id
 from quantdsl.domain.model.contract_valuation import ContractValuation
 from quantdsl.domain.model.market_simulation import MarketSimulation
-from quantdsl.domain.model.simulated_price import make_simulated_price_id
+from quantdsl.domain.model.simulated_price import make_simulated_price_id, SimulatedPrice
 from quantdsl.domain.model.simulated_price_requirements import SimulatedPriceRequirements
 from quantdsl.exceptions import TimeoutError, InterruptSignalReceived
 from quantdsl.priceprocess.base import datetime_from_date, get_duration_years
+from quantdsl.semantics import discount
 
 
 class Results(object):
@@ -93,10 +95,11 @@ def calc_print(source_code, observation_date=None, interest_rate=0, path_count=2
     return results
 
 
-def calc(source_code, observation_date=None, interest_rate=0, path_count=20000, perturbation_factor=0.01,
-         price_process=None, periodisation=None, dsl_classes=None,
+def calc(source_code, observation_date=None, interest_rate=0, path_count=20000,
+         perturbation_factor=0.01, price_process=None, periodisation=None, dsl_classes=None,
          max_dependency_graph_size=DEFAULT_MAX_DEPENDENCY_GRAPH_SIZE,
          timeout=None, verbose=False):
+
     cmd = Calculate(
         source_code=source_code,
         observation_date=observation_date,
@@ -313,24 +316,39 @@ def read_results(self, app, evaluation, market_simulation):
                 if len(perturbed_name_split) > 3:
                     day = int(perturbed_name_split[3])
                     price_date = datetime.date(year, month, day)
-                else:
-                    price_date = datetime.date(year, month, 1)
-                simulated_price_id = make_simulated_price_id(market_simulation.id, commodity_name, price_date,
-                                                             price_date)
-                try:
+                    simulated_price_id = make_simulated_price_id(
+                        market_simulation.id, commodity_name, price_date, price_date
+                    )
                     simulated_price = app.simulated_price_repo[simulated_price_id]
-                except KeyError as e:
-                    raise Exception("Simulated price for date {} is unavailable".format(price_date, e))
+                    simulated_price_value = simulated_price.value
+                else:
+                    sum_simulated_prices = 0
+                    count_simulated_prices = 0
+                    for i in range(1, 32):
+                        try:
+                            price_date = datetime.date(year, month, i)
+                        except ValueError:
+                            continue
+                        else:
+                            simulated_price_id = make_simulated_price_id(
+                                market_simulation.id, commodity_name, price_date, price_date
+                            )
+                            try:
+                                simulated_price = app.simulated_price_repo[simulated_price_id]
+                            except KeyError:
+                                pass
+                            else:
+                                sum_simulated_prices += simulated_price.value
+                                count_simulated_prices += 1
+                    assert count_simulated_prices, "Can't find any simulated prices for {}-{}".format(year, month)
+                    simulated_price_value = sum_simulated_prices / count_simulated_prices
 
                 dy = perturbed_value - perturbed_value_negative
-                price = simulated_price.value
+                price = simulated_price_value
                 dx = 2 * market_simulation.perturbation_factor * price
                 contract_delta = dy / dx
-                years = get_duration_years(market_simulation.observation_date, price_date)
                 # Todo: Refactor this w.r.t the discount() method of DslExpression.
-                interest_rate = market_simulation.interest_rate / 100.0
-                discount_rate = math.exp(-interest_rate * years)
-                hedge_units = - contract_delta / discount_rate
+                hedge_units = - discount(contract_delta, price_date, market_simulation.observation_date, market_simulation.interest_rate)
                 cash_in = contract_delta * price
                 periods.append({
                     'commodity': perturbed_name,
@@ -554,6 +572,11 @@ def plot_periods(periods, title, periodisation, interest_rate, path_count, pertu
     else:
         return
 
+    NUM_STD_DEVS = 2
+    OUTER_COLOUR = 'y'
+    MID_COLOUR = 'g'
+    INNER_COLOUR = 'b'
+    MEAN_COLOUR = '0.1'
     for i, name in enumerate(names):
 
         _periods = [p for p in periods if p['commodity'].startswith(name)]
@@ -567,8 +590,8 @@ def plot_periods(periods, title, periodisation, interest_rate, path_count, pertu
         # prices_3 = [p['price'][2] for p in _periods]
         # prices_4 = [p['price'][3] for p in _periods]
         prices_std = [p['price'].std() for p in _periods]
-        prices_plus = list(numpy.array(prices_mean) + 2 * numpy.array(prices_std))
-        prices_minus = list(numpy.array(prices_mean) - 2 * numpy.array(prices_std))
+        prices_plus = list(numpy.array(prices_mean) + NUM_STD_DEVS * numpy.array(prices_std))
+        prices_minus = list(numpy.array(prices_mean) - NUM_STD_DEVS * numpy.array(prices_std))
 
         price_plot.set_title('Prices - {}'.format(name))
         price_plot.plot(
@@ -592,22 +615,37 @@ def plot_periods(periods, title, periodisation, interest_rate, path_count, pertu
                 pos += cum_pos[-1]
             cum_pos.append(pos)
         cum_pos_mean = [p.mean() for p in cum_pos]
+        cum_pos_p5 = [nanpercentile(p, 5) for p in cum_pos]
+        cum_pos_p10 = [nanpercentile(p, 10) for p in cum_pos]
+        cum_pos_p25 = [nanpercentile(p, 25) for p in cum_pos]
+        cum_pos_p75 = [nanpercentile(p, 75) for p in cum_pos]
+        cum_pos_p90 = [nanpercentile(p, 90) for p in cum_pos]
+        cum_pos_p95 = [nanpercentile(p, 95) for p in cum_pos]
+
         cum_pos_std = [p.std() for p in cum_pos]
         cum_pos_stderr = [p / math.sqrt(path_count) for p in cum_pos_std]
-        cum_pos_std_plus = list(numpy.array(cum_pos_mean) + 1 * numpy.array(cum_pos_std))
-        cum_pos_std_minus = list(numpy.array(cum_pos_mean) - 1 * numpy.array(cum_pos_std))
-        cum_pos_stderr_plus = list(numpy.array(cum_pos_mean) + 3 * numpy.array(cum_pos_stderr))
-        cum_pos_stderr_minus = list(numpy.array(cum_pos_mean) - 3 * numpy.array(cum_pos_stderr))
+        cum_pos_std_offset = NUM_STD_DEVS * numpy.array(cum_pos_std)
+        cum_pos_std_plus = list(numpy.array(cum_pos_mean) + cum_pos_std_offset)
+        cum_pos_std_minus = list(numpy.array(cum_pos_mean) - cum_pos_std_offset)
+        com_pos_stderr_offset = NUM_STD_DEVS * numpy.array(cum_pos_stderr)
+        cum_pos_stderr_plus = list(numpy.array(cum_pos_mean) + com_pos_stderr_offset)
+        cum_pos_stderr_minus = list(numpy.array(cum_pos_mean) - com_pos_stderr_offset)
 
         pos_plot = subplots[len(names) + i]
         pos_plot.set_title('Position - {}'.format(name))
 
         pos_plot.plot(
-            dates, cum_pos_std_plus, '0.85',
-            dates, cum_pos_std_minus, '0.85',
-            dates, cum_pos_stderr_plus, '0.5',
-            dates, cum_pos_stderr_minus, '0.5',
-            dates, cum_pos_mean, '0.1',
+            # dates, cum_pos_std_plus, '0.85',
+            # dates, cum_pos_std_minus, '0.85',
+            dates, cum_pos_p5, OUTER_COLOUR,
+            dates, cum_pos_p10, MID_COLOUR,
+            dates, cum_pos_p25, INNER_COLOUR,
+            dates, cum_pos_p75, INNER_COLOUR,
+            dates, cum_pos_p90, MID_COLOUR,
+            dates, cum_pos_p95, OUTER_COLOUR,
+            # dates, cum_pos_stderr_plus, '0.5',
+            # dates, cum_pos_stderr_minus, '0.5',
+            dates, cum_pos_mean, MEAN_COLOUR,
         )
 
         ymin = min(0, min(cum_pos_std_minus)) - 1
@@ -633,20 +671,36 @@ def plot_periods(periods, title, periodisation, interest_rate, path_count, pertu
             cash_in += cum_cash_in[-1]
         cum_cash_in.append(cash_in)
 
+    cum_cash_p5 = [nanpercentile(p, 5) for p in cum_cash_in]
+    cum_cash_p10 = [nanpercentile(p, 10) for p in cum_cash_in]
+    cum_cash_p25 = [nanpercentile(p, 25) for p in cum_cash_in]
+    cum_cash_p75 = [nanpercentile(p, 75) for p in cum_cash_in]
+    cum_cash_p90 = [nanpercentile(p, 90) for p in cum_cash_in]
+    cum_cash_p95 = [nanpercentile(p, 95) for p in cum_cash_in]
     cum_cash_mean = [p.mean() for p in cum_cash_in]
     cum_cash_std = [p.std() for p in cum_cash_in]
     cum_cash_stderr = [p / math.sqrt(path_count) for p in cum_cash_std]
-    cum_cash_std_plus = list(numpy.array(cum_cash_mean) + 1 * numpy.array(cum_cash_std))
-    cum_cash_std_minus = list(numpy.array(cum_cash_mean) - 1 * numpy.array(cum_cash_std))
-    cum_cash_stderr_plus = list(numpy.array(cum_cash_mean) + 3 * numpy.array(cum_cash_stderr))
-    cum_cash_stderr_minus = list(numpy.array(cum_cash_mean) - 3 * numpy.array(cum_cash_stderr))
+
+    cum_cash_std_offset = NUM_STD_DEVS * numpy.array(cum_cash_std)
+    cum_cash_std_plus = list(numpy.array(cum_cash_mean) + cum_cash_std_offset)
+    cum_cash_std_minus = list(numpy.array(cum_cash_mean) - cum_cash_std_offset)
+
+    cum_cash_stderr_offset = NUM_STD_DEVS * numpy.array(cum_cash_stderr)
+    cum_cash_stderr_plus = list(numpy.array(cum_cash_mean) + cum_cash_stderr_offset)
+    cum_cash_stderr_minus = list(numpy.array(cum_cash_mean) - cum_cash_stderr_offset)
 
     profit_plot.plot(
-        dates, cum_cash_std_plus, '0.8',
-        dates, cum_cash_std_minus, '0.8',
-        dates, cum_cash_stderr_plus, '0.5',
-        dates, cum_cash_stderr_minus, '0.5',
-        dates, cum_cash_mean, '0.1'
+        # dates, cum_cash_std_plus, '0.8',
+        # dates, cum_cash_std_minus, '0.8',
+        dates, cum_cash_p5, OUTER_COLOUR,
+        dates, cum_cash_p10, MID_COLOUR,
+        dates, cum_cash_p25, INNER_COLOUR,
+        dates, cum_cash_p75, INNER_COLOUR,
+        dates, cum_cash_p90, MID_COLOUR,
+        dates, cum_cash_p95, OUTER_COLOUR,
+        # dates, cum_cash_stderr_plus, '0.5',
+        # dates, cum_cash_stderr_minus, '0.5',
+        dates, cum_cash_mean, MEAN_COLOUR
     )
 
     f.autofmt_xdate(rotation=60)

quantdsl/semantics.py

@@ -213,6 +213,16 @@ def identify_perturbation_dependencies(self, dependencies, **kwds):
                 dsl_arg.identify_perturbation_dependencies(dependencies, **kwds)
 
 
+def discount(value, start_date, end_date, interest_rate):
+    r = interest_rate / 100
+    T = get_duration_years(start_date, end_date)
+    # Assumes continuous compounding rate.
+    discount_factor = math.exp(- r * T)
+    return value * discount_factor
+    # Not annualised equivalent rate.
+    # return value / (1 + r) ** T
+
+
 class DslExpression(DslObject):
     relative_cost = 1
 
@@ -222,15 +232,6 @@ def evaluate(self, **kwds):
         Returns the value of the expression.
         """
 
-    def discount(self, value, start_date, end_date, **kwds):
-        r = float(kwds['interest_rate']) / 100
-        T = get_duration_years(start_date, end_date)
-        # Assumes continuous compounding rate.
-        discount_factor = math.exp(- r * T)
-        return value * discount_factor
-        # Not annualised equivalent rate.
-        # return value / (1 + r) ** T
-
     def cost_element(self):
         return self.relative_cost
 
@@ -1392,7 +1393,8 @@ def evaluate(self, **kwds):
         included_value = self._args[1].evaluate(**kwds)
 
         present_time = self.get_present_time(kwds)
-        discounted_value = self.discount(included_value, present_time, self.get_date(**kwds), **kwds)
+        interest_rate = kwds['interest_rate']
+        discounted_value = discount(included_value, present_time, self.get_date(**kwds), interest_rate)
         return discounted_value
 
 
@@ -1463,6 +1465,7 @@ class On(Fixing):
 
 
 class Wait(Fixing):
+# class Wait(Settlement, Fixing):
     """
     A fixing with discounting of the resulting value from date arg to present_time.
 
@@ -1471,9 +1474,11 @@ class Wait(Fixing):
     relative_cost = 10
 
     def evaluate(self, **kwds):
+        # return super(Wait, self).evaluate(**kwds)
         value = super(Wait, self).evaluate(**kwds)
         present_time = self.get_present_time(kwds)
-        return self.discount(value, present_time, self.get_date(**kwds), **kwds)
+        interest_rate = kwds['interest_rate']
+        return discount(value, present_time, self.get_date(**kwds), interest_rate)
 
 
 class Choice(StochasticObject, DslExpression):