Merge pull request #5 from EconomicSL/master

Non-prpoportional

distribution_wrapper_test.py

@@ -0,0 +1,39 @@
+import scipy.stats
+import numpy as np
+import matplotlib.pyplot as plt
+from distributiontruncated import TruncatedDistWrapper
+from distributionreinsurance import ReinsuranceDistWrapper
+import pdb
+
+non_truncated_dist = scipy.stats.pareto(b=2, loc=0, scale=0.5)
+truncated_dist = TruncatedDistWrapper(lower_bound=0.6, upper_bound=1., dist=non_truncated_dist)
+reinsurance_dist = ReinsuranceDistWrapper(lower_bound=0.85, upper_bound=0.95, dist=truncated_dist)
+
+x1 = np.linspace(non_truncated_dist.ppf(0.01), non_truncated_dist.ppf(0.99), 100)
+x2 = np.linspace(truncated_dist.ppf(0.01), truncated_dist.ppf(1.), 100)
+x3 = np.linspace(reinsurance_dist.ppf(0.01), reinsurance_dist.ppf(1.), 100)
+x_val_1 = reinsurance_dist.lower_bound
+x_val_2 = truncated_dist.upper_bound - (reinsurance_dist.upper_bound - reinsurance_dist.lower_bound)
+x_val_3 = reinsurance_dist.upper_bound
+x_val_4 = truncated_dist.upper_bound
+
+fig, ax = plt.subplots(1, 1)
+ax.plot(x1, non_truncated_dist.pdf(x1), 'k-', lw=2, label='non-truncated pdf')
+ax.plot(x1, non_truncated_dist.cdf(x1), 'g-', lw=2, label='non-truncated cdf')
+ax.plot(x2, truncated_dist.pdf(x2), 'r-', lw=2, label='truncated pdf')
+ax.plot(x2, truncated_dist.cdf(x2), 'm-', lw=2, label='truncated cdf')
+ax.plot(x3, reinsurance_dist.pdf(x3), 'b-', lw=2, label='reinsurance pdf')
+ax.plot(x3, reinsurance_dist.cdf(x3), 'c-', lw=2, label='reinsurance cdf')
+ax.set_xlim(0.45, 1.25)
+ax.set_ylim(0, 5)
+ax.arrow(x_val_1, reinsurance_dist.cdf(x_val_1), x_val_3 - x_val_1, truncated_dist.cdf(x_val_3) - reinsurance_dist.cdf(x_val_1), head_width=0, head_length=0, fc='m', ec='m', ls=':')
+ax.arrow(x_val_2, reinsurance_dist.cdf(x_val_2), x_val_4 - x_val_2, truncated_dist.cdf(x_val_4) - reinsurance_dist.cdf(x_val_2), head_width=0, head_length=0, fc='m', ec='m', ls=':')
+ax.arrow(x_val_1, reinsurance_dist.pdf(x_val_1+0.00001), x_val_3 - x_val_1, truncated_dist.pdf(x_val_3) - reinsurance_dist.pdf(x_val_1+0.00001), head_width=0, head_length=0, fc='r', ec='r', ls=':')
+ax.arrow(x_val_2, reinsurance_dist.pdf(x_val_2), x_val_4 - x_val_2, truncated_dist.pdf(x_val_4) - reinsurance_dist.pdf(x_val_2), head_width=0, head_length=0, fc='r', ec='r', ls=':')
+sample = reinsurance_dist.rvs(size=100000)
+#sample = sample[sample < scipy.percentile(sample, 90)]
+ax.hist(sample, normed=True, histtype='stepfilled', alpha=0.4)
+ax.legend(loc='best', frameon=False)
+plt.show()
+
+pdb.set_trace()

distributionreinsurance.py

@@ -0,0 +1,85 @@
+import scipy.stats
+import numpy as np
+import matplotlib.pyplot as plt
+from math import ceil
+import scipy
+import pdb
+
+class ReinsuranceDistWrapper():
+    def __init__(self, dist, lower_bound=None, upper_bound=None):
+        assert lower_bound is not None or upper_bound is not None
+        self.dist = dist
+        self.lower_bound = lower_bound
+        self.upper_bound = upper_bound
+        if lower_bound is None:
+            self.lower_bound = -np.inf
+        elif upper_bound is None:
+            self.upper_bound = np.inf
+        assert self.upper_bound > self.lower_bound
+        self.redistributed_share = dist.cdf(upper_bound) - dist.cdf(lower_bound)
+
+
+    def pdf(self, x):
+        x = np.array(x, ndmin=1)
+        r = map(lambda Y: self.dist.pdf(Y) if Y < self.lower_bound \
+                            else np.inf if Y==self.lower_bound \
+                            else self.dist.pdf(Y + self.upper_bound - self.lower_bound), x)
+        r = np.array(list(r))
+        if len(r.flatten()) == 1:
+            r = float(r)
+        return r
+
+    def cdf(self, x):
+        x = np.array(x, ndmin=1)
+        r = map(lambda Y: self.dist.cdf(Y) if Y < self.lower_bound \
+                        else self.dist.cdf(Y + self.upper_bound - self.lower_bound), x)
+        r = np.array(list(r))
+        if len(r.flatten()) == 1:
+            r = float(r)
+        return r
+
+    def ppf(self, x):
+        x = np.array(x, ndmin=1)
+        assert (x >= 0).all() and (x <= 1).all()
+        r = map(lambda Y: self.dist.ppf(Y) if Y <= self.dist.cdf(self.lower_bound) \
+                        else self.dist.ppf(self.dist.cdf(self.lower_bound)) if Y <= self.dist.cdf(self.upper_bound) \
+                        else self.dist.ppf(Y) - self.upper_bound + self.lower_bound, x)
+        r = np.array(list(r))
+        if len(r.flatten()) == 1:
+            r = float(r)
+        return r
+
+    def rvs(self, size=1):
+        sample = self.dist.rvs(size=size)
+        sample1 = sample[sample<=self.lower_bound]
+        sample2 = sample[sample>self.lower_bound]
+        sample3 = sample2[sample2>=self.upper_bound]
+        sample2 = sample2[sample2<self.upper_bound]
+
+        sample2 = np.ones(len(sample2)) * self.lower_bound
+        sample3 = sample3 -self.upper_bound + self.lower_bound
+
+        sample = np.append(np.append(sample1,sample2),sample3)
+        return sample[:size]
+
+
+if __name__ == "__main__":
+    non_truncated = scipy.stats.pareto(b=2, loc=0, scale=0.5)
+    #truncated = ReinsuranceDistWrapper(lower_bound=0, upper_bound=1, dist=non_truncated)
+    truncated = ReinsuranceDistWrapper(lower_bound=0.9, upper_bound=1.1, dist=non_truncated)
+
+    x = np.linspace(non_truncated.ppf(0.01), non_truncated.ppf(0.99), 100)
+    x2 = np.linspace(truncated.ppf(0.01), truncated.ppf(0.99), 100)
+
+    fig, ax = plt.subplots(1, 1)
+    ax.plot(x, non_truncated.pdf(x), 'k-', lw=2, label='non-truncated pdf')
+    ax.plot(x2, truncated.pdf(x2), 'r-', lw=2, label='truncated pdf')
+    ax.plot(x, non_truncated.cdf(x), 'b-', lw=2, label='non-truncated cdf')
+    ax.plot(x2, truncated.cdf(x2), 'g-', lw=2, label='truncated cdf')
+    sample = truncated.rvs(size=1000)
+    sample = sample[sample < scipy.percentile(sample, 90)]
+    ax.hist(sample, normed=True, histtype='stepfilled', alpha=0.4)
+    ax.legend(loc='best', frameon=False)
+    plt.show()
+
+    #pdb.set_trace()

distributiontruncated.py

@@ -0,0 +1,68 @@
+import scipy.stats
+import numpy as np
+import matplotlib.pyplot as plt
+from math import ceil
+import scipy.integrate
+
+class TruncatedDistWrapper():
+    def __init__(self, dist, lower_bound=0, upper_bound=1):
+        self.dist = dist
+        self.normalizing_factor = dist.cdf(upper_bound) - dist.cdf(lower_bound)
+        self.lower_bound = lower_bound
+        self.upper_bound = upper_bound
+        assert self.upper_bound > self.lower_bound
+
+    def pdf(self, x):
+        x = np.array(x, ndmin=1)
+        r = map(lambda Y: self.dist.pdf(Y) / self.normalizing_factor \
+                            if (Y >= self.lower_bound and Y <= self.upper_bound) else 0, x)
+        r = np.array(list(r))
+        if len(r.flatten()) == 1:
+            r = float(r)
+        return r
+
+    def cdf(self, x):
+        x = np.array(x, ndmin=1)
+        r = map(lambda Y: 0 if Y < self.lower_bound else 1 if Y > self.upper_bound \
+                            else (self.dist.cdf(Y) - self.dist.cdf(self.lower_bound))/ self.normalizing_factor, x)
+        r = np.array(list(r))
+        if len(r.flatten()) == 1:
+            r = float(r)
+        return r
+
+    def ppf(self, x):
+        x = np.array(x, ndmin=1)
+        assert (x >= 0).all() and (x <= 1).all()
+        return self.dist.ppf(x * self.normalizing_factor + self.dist.cdf(self.lower_bound))
+
+    def rvs(self, size=1):
+        init_sample_size = int(ceil(size / self.normalizing_factor * 1.1))
+        sample = self.dist.rvs(size=init_sample_size)
+        sample = sample[sample>=self.lower_bound]
+        sample = sample[sample<=self.upper_bound]
+        while len(sample) < size:
+            sample = np.append(sample, self.rvs(size - len(sample)))
+        return sample[:size]    
+
+    def mean(self):
+        mean_estimate, mean_error = scipy.integrate.quad(lambda Y: Y*self.pdf(Y), self.lower_bound, self.upper_bound)
+        return mean_estimate
+
+if __name__ == "__main__":
+    non_truncated = scipy.stats.pareto(b=2, loc=0, scale=0.5)
+    truncated = TruncatedDistWrapper(lower_bound=0.55, upper_bound=1., dist=non_truncated)
+
+    x = np.linspace(non_truncated.ppf(0.01), non_truncated.ppf(0.99), 100)
+    x2 = np.linspace(truncated.ppf(0.01), truncated.ppf(0.99), 100)
+
+    fig, ax = plt.subplots(1, 1)
+    ax.plot(x, non_truncated.pdf(x), 'k-', lw=2, label='non-truncated pdf')
+    ax.plot(x2, truncated.pdf(x2), 'r-', lw=2, label='truncated pdf')
+    ax.plot(x, non_truncated.cdf(x), 'b-', lw=2, label='non-truncated cdf')
+    ax.plot(x2, truncated.cdf(x2), 'g-', lw=2, label='truncated cdf')
+    sample = truncated.rvs(size=1000)
+    ax.hist(sample, normed=True, histtype='stepfilled', alpha=0.4)
+    ax.legend(loc='best', frameon=False)
+    plt.show()
+
+    print(truncated.mean())

insurancecontract.py

@@ -2,7 +2,7 @@
 import sys, pdb
 
 class InsuranceContract():
-    def __init__(self, insurer, properties, time, premium, runtime, payment_period, deductible=0, excess=None, reinsurance=0):
+    def __init__(self, insurer, properties, time, premium, runtime, payment_period, insurancetype="proportional", deductible=0, excess=None, reinsurance=0):
         """Constructor method.
                Accepts arguments
                     insurer: Type InsuranceFirm. 
@@ -12,6 +12,7 @@ def __init__(self, insurer, properties, time, premium, runtime, payment_period,
                     runtime: Type integer.
                     payment_period: Type integer.
                 optional:
+                    insurancetype: Type string. The type of this contract, especially "proportional" vs "excess_of_loss"
                     deductible: Type float (or int)
                     excess: Type float (or int or None)
                     reinsurance: Type float (or int). The value that is being reinsured.
@@ -27,8 +28,10 @@ def __init__(self, insurer, properties, time, premium, runtime, payment_period,
         self.property_holder = properties["owner"]
         self.value = properties["value"]
         self.contract = properties.get("contract")  # will assign None if key does not exist
+        self.insurancetype = properties.get("insurancetype") if insurancetype is None else insurancetype 
         self.properties = properties
         self.runtime = runtime
+        self.starttime = time
         self.expiration = runtime + time
         self.terminating = False
 
@@ -38,7 +41,7 @@ def __init__(self, insurer, properties, time, premium, runtime, payment_period,
         #    self.deductible = deductible if deductible is not None else 0
         self.deductible = deductible
 
-        self.excess = excess if excess is not None else self.value
+        self.excess = excess if excess is not None else 1.
 
         self.reinsurance = reinsurance
         self.reinsurer = None
@@ -47,8 +50,11 @@ def __init__(self, insurer, properties, time, premium, runtime, payment_period,
         #self.is_reinsurancecontract = False
 
         # setup payment schedule
+        #total_premium = premium * (self.excess - self.deductible)   # TODO: excess and deductible should not be considered linearily in premium computation; this should be shifted to the (re)insurer who supplies the premium as argument to the contract's constructor method
+        total_premium = premium * self.value 
+        self.periodized_premium = total_premium / self.runtime
         self.payment_times = [time + i for i in range(runtime) if i % payment_period == 0]
-        self.payment_values = premium * (np.ones(len(self.payment_times)) / len(self.payment_times))
+        self.payment_values = total_premium * (np.ones(len(self.payment_times)) / len(self.payment_times))
 
         ## Create obligation for premium payment
         #self.property_holder.receive_obligation(premium * (self.excess - self.deductible), self.insurer, time)
@@ -62,7 +68,7 @@ def check_payment_due(self, time):
         if len(self.payment_times) > 0 and time >= self.payment_times[0]:
             # Create obligation for premium payment
             #self.property_holder.receive_obligation(premium * (self.excess - self.deductible), self.insurer, time)
-            self.property_holder.receive_obligation(self.payment_values[0] * (self.excess - self.deductible), self.insurer, time)
+            self.property_holder.receive_obligation(self.payment_values[0], self.insurer, time)
 
             # Remove current payment from payment schedule
             self.payment_times = self.payment_times[1:]

insurancefirm.py

@@ -31,14 +31,20 @@ def init(self, simulation_parameters, agent_parameters):
         self.acceptance_threshold_friction = agent_parameters['acceptance_threshold_friction']  # 0.9 #1.0 to switch off
         self.interest_rate = agent_parameters["interest_rate"]
         self.reinsurance_limit = agent_parameters["reinsurance_limit"]
+        self.simulation_no_risk_categories = simulation_parameters["no_categories"]
+        self.simulation_reinsurance_type = simulation_parameters["simulation_reinsurance_type"]
+        self.categories_reinsured = [False for i in range(self.simulation_no_risk_categories)]
+        if self.simulation_reinsurance_type == 'non-proportional':
+            self.np_reinsurance_deductible = simulation_parameters["default_non-proportional_reinsurance_deductible"]
+            self.np_reinsurance_excess = simulation_parameters["default_non-proportional_reinsurance_excess"]
         self.obligations = []
         self.underwritten_contracts = []
         #self.reinsurance_contracts = []
         self.operational = True
         self.is_insurer = True
         self.is_reinsurer = False
 
-    def iterate(self, time):
+    def iterate(self, time):        # TODO: split function so that only the sequence of events remains here and everything else is in separate methods
         """obtain investments yield"""
         self.obtain_yield(time)
 
@@ -71,12 +77,24 @@ def iterate(self, time):
             except:
                 print("Something wrong; agent {0:d} receives too few new contracts {1:d} <= {2:d}".format(self.id, contracts_offered, 2*contracts_dissolved))
             #print(self.id, " has ", len(self.underwritten_contracts), " & receives ", contracts_offered, " & lost ", contracts_dissolved)
+
+            new_nonproportional_risks = [risk for risk in new_risks if risk.get("insurancetype")=='excess-of-loss']
+            new_risks = [risk for risk in new_risks if risk.get("insurancetype") in ['proportional', None]]
 
-
-            """make underwriting decisions, category-wise"""
             underwritten_risks = [{"excess": contract.value, "category": contract.category, \
                             "risk_factor": contract.risk_factor, "deductible": contract.deductible, \
+                            "excess": contract.excess, "insurancetype": contract.insurancetype, \
                             "runtime": contract.runtime} for contract in self.underwritten_contracts if contract.reinsurance_share != 1.0]
+
+            """deal with non-proportional risks first as they must evaluate each request separatly, then with proportional ones"""
+            for risk in new_nonproportional_risks:
+                #accept = self.riskmodel.evaluate(underwritten_risks, self.cash, risk)       # TODO: change riskmodel.evaluate() to accept new risk to be evaluated and to account for existing non-proportional risks correctly
+                #if accept:
+                #    contract = ReinsuranceContract(...)
+                #    self.underwritten_contracts.append(contract)
+                pass    # TODO: write this nonproportional risk acceptance decision section based on commented code in the lines above this
+
+            """make underwriting decisions, category-wise"""
             # TODO: Enable reinsurance shares other tan 0.0 and 1.0
             expected_profit, acceptable_by_category = self.riskmodel.evaluate(underwritten_risks, self.cash)
 
@@ -171,9 +189,41 @@ def receive(self, amount):
     def obtain_yield(self, time):
         amount = self.cash * self.interest_rate
         self.simulation.receive_obligation(amount, self, time)
+
+    def ask_reinsurance(self):
+        if self.simulation_reinsurance_type == 'proportional':
+            self.ask_reinsurance_proportional()
+        elif self.simulation_reinsurance_type == 'non-proportional':
+            self.ask_reinsurance_non_proportional()
+        else:
+            assert False, "Undefined reinsurance type"
+
+    def ask_reinsurance_non_proportional(self):
+        for categ_id in range(len(self.simulation_no_risk_categories)):
+            # with probability 5% if not reinsured ...      # TODO: find a more generic way to decide whether to request reinsurance for category in this period
+            if (not self.categories_reinsured[categ_id]) and np.random() < 0.05:
+                total_value = 0
+                avg_risk_factor = 0
+                number_risks = 0
+                periodized_total_premium = 0
+                for contract in self.underwritten_contracts:
+                    if contract.category == categ_id:
+                        total_value += contract.value
+                        avg_risk_factor += contract.risk_factor
+                        number_risks += 1
+                        periodized_total_premium += contract.periodized_premium
+                avg_risk_factor /= number_risks
+                risk = {"value": total_value, "category": categ_id, "owner": self,
+                            #"identifier": uuid.uuid1(),
+                            "insurancetype": 'excess-of-loss', "number_risks": number_risks, 
+                            "deductible": self.np_reinsurance_deductible, "excess": np_reinsurance_excess,
+                            "periodized_total_premium": periodized_total_premium, "runtime": 12,
+                            "expiration": time + 12, "risk_factor": avg_risk_factor}
+
+                self.simulation.append_reinrisks(risk)
 
     @nb.jit
-    def ask_reinsurance(self):
+    def ask_reinsurance_proportional(self):
         nonreinsured = []
         for contract in self.underwritten_contracts:
             if contract.reincontract == None:

insurancesimulation.py

@@ -2,6 +2,7 @@
 from insurancefirm import InsuranceFirm
 from riskmodel import RiskModel
 from reinsurancefirm import ReinsuranceFirm
+from distributiontruncated import TruncatedDistWrapper
 import numpy as np
 import scipy.stats
 import math
@@ -32,7 +33,14 @@ def __init__(self, override_no_riskmodels, replic_ID, simulation_parameters):
         self.simulation_parameters = simulation_parameters
 
         # unpack parameters, set up environment (distributions etc.)
-        self.damage_distribution = scipy.stats.uniform(loc=0, scale=1)
+
+        # damage distribution
+        # TODO: control damage distribution via parameters, not directly
+        #self.damage_distribution = scipy.stats.uniform(loc=0, scale=1)
+        non_truncated = scipy.stats.pareto(b=2, loc=0, scale=0.25)
+        self.damage_distribution = TruncatedDistWrapper(lower_bound=0.25, upper_bound=1., dist=non_truncated)
+
+        # remaining parameters
         self.cat_separation_distribution = scipy.stats.expon(0, simulation_parameters["event_time_mean_separation"])
         self.risk_factor_lower_bound = simulation_parameters["risk_factor_lower_bound"]
         self.risk_factor_spread = simulation_parameters["risk_factor_upper_bound"] - simulation_parameters["risk_factor_lower_bound"]

isleconfig.py

@@ -21,6 +21,9 @@
                        "risk_factor_upper_bound": 0.6,
                        "initial_acceptance_threshold": 0.5,
                        "acceptance_threshold_friction": 0.9,
+                       "simulation_reinsurance_type": 'proportional',
+                       "default_non-proportional_reinsurance_deductible": 0.75,
+                       "default_non-proportional_reinsurance_excess": 1.0,
                        "initial_agent_cash": 10000,
                        "initial_reinagent_cash": 50000,
                        "interest_rate": 0,