import pandas as pd
import numpy as np
import chainladder as cl
import seaborn as sns
sns.set_style('whitegrid')
%matplotlib inline
print('chainladder:' + cl.__version__)
print('pandas:' + pd.__version__)

"""
Explanation: Stochastic methods
Getting started
All exercises rely on chainladder v0.5.5 and later.  There have also been breaking changes with pandas 1.0 and if you are using an earlier version, date slicing may behave differently.
End of explanation
"""


tri = cl.load_sample('clrd').groupby('LOB').sum().loc['wkcomp', ['CumPaidLoss', 'EarnedPremNet']]
cl.Chainladder().fit(tri['CumPaidLoss']).ultimate_ == \
cl.MackChainladder().fit(tri['CumPaidLoss']).ultimate_

"""
Explanation: MackChainladder Intro
Like the basic Chainladder method, the MackChainladder is entirely specified by its development pattern selections.  In fact, it is the basic Chainladder with a few extra features.  Let's explore this a bit more with the Workers' Compensation industry triangle.
End of explanation
"""


mack = cl.MackChainladder().fit(tri['CumPaidLoss'])

"""
Explanation: Let's create a Mack Model.
End of explanation
"""


tri_first_lags = tri[tri.development<=24][tri.origin<'1997']['CumPaidLoss']
tri_first_lags

"""
Explanation: MackChainladder has the following additional fitted features that the deterministic Chainladder does not.

full_std_err_:  The full standard error
total_process_risk_: The total process error
total_parameter_risk_: The total parameter error
mack_std_err_: The total prediction error by origin period
total_mack_std_err_: The total prediction error across all origin periods

Notice these are all measures of uncertainty, but where do they come from?  Let's start by examining the link_ratios underlying the triangle.
End of explanation
"""


tri_first_lags.link_ratio.to_frame().mean().values[0]

"""
Explanation: A simple average link-ratio can be directly computed as follows:
End of explanation
"""


cl.Development(average='simple').fit(tri['CumPaidLoss']).ldf_.to_frame().values[0, 0]

"""
Explanation: Verifying that this ties to our Development object:
End of explanation
"""


import statsmodels.api as sm
import numpy as np
y = tri_first_lags.to_frame().values[:, 1]
X = tri_first_lags.to_frame().values[:, 0]

model = sm.WLS(y, X, weights=(1/X)**2)
results = model.fit()
results.summary()

"""
Explanation: A linear regression framework
Mack noticed that this estimate for an LDF is really just a linear regression fit.  For the case of the simple average, it is a weighted regression where the weight is set to $\left (\frac{1}{X}  \right )^{2}$.
Take a look at the fitted coefficient in the next cell and verify that it ties to the direct calculations above.
With the regression framework in hand, we get much more information about our LDF estimate than just the coefficient.
End of explanation
"""


print('Does this work for simple?')
print(round(cl.Development(average='simple').fit(tri_first_lags).ldf_.to_frame().values[0, 0], 8) == \
      round(sm.WLS(y, X, weights=(1/X)**2).fit().params[0],8))
print('Does this work for volume-weighted average?')
print(round(cl.Development(average='volume').fit(tri_first_lags).ldf_.to_frame().values[0, 0], 8) == \
      round(sm.WLS(y, X, weights=(1/X)).fit().params[0],8))
print('Does this work for regression average?')
print(round(cl.Development(average='regression').fit(tri_first_lags).ldf_.to_frame().values[0, 0], 8) == \
      round(sm.OLS(y, X).fit().params[0],8))

"""
Explanation: By toggling the weights of our regression, we can handle the most common types of averaging used in picking loss development factors.
End of explanation
"""


dev = cl.Development(average='simple').fit(tri['CumPaidLoss'])

dev.std_err_

dev.sigma_

"""
Explanation: This regression framework is what the Development estimator uses to set development patterns.  Although we discard the information in deterministic approaches, Development has two useful statistics for estimating reserve variability, both of which come from the regression framework.  The stastics are std_err_ and sigma_ and they are used by the MackChainladder estimator to determine the prediction error of our reserves.
End of explanation
"""


print('Does this work for dropping observations?')
print(round(cl.Development(average='volume', drop_valuation='1988') \
              .fit(tri['CumPaidLoss']).std_err_.to_frame().values[0, 0], 8) == \
      round(sm.WLS(y[1:], X[1:], weights=(1/X[1:])).fit().bse[0],8))

"""
Explanation: Since the regression framework is weighted, we can easily turn on/off any observation we want using the dropping capabilities of the Development estimator.  Dropping link ratios not only affects the ldf_ and cdf_, but also the std_err_ and sigma of the regression.
Here we eliminate the 1988 valuation from our triangle, which is identical to eliminating the first observation from our 12-24 regression fit.
End of explanation
"""


mack.parameter_risk_

"""
Explanation: With sigma_ and std_err_ in hand, Mack goes on to develop recursive formulas to estimate parameter_risk_ and process_risk_.
End of explanation
"""


print('Parameter risk and process risk are independent?')
print(round(mack.mack_std_err_**2, 4) == round(mack.parameter_risk_**2 + mack.process_risk_**2, 4))

"""
Explanation: Assumption of Independence
The Mack model makes a lot of assumptions about independence (i.e. covariance between random processes is 0).  This means many of the Variance estimates in the MackChainladder model follow the form of $Var(A+B) = Var(A)+Var(B)$.
Notice the square of mack_std_err_ is simply the sum of the sqaures of parameter_risk_ and process_risk_.
End of explanation
"""


print('Total Parameter and process risk across origin periods is independent?')
print(round(mack.total_process_risk_**2, 4) == round((mack.process_risk_**2).sum('origin'), 4))

"""
Explanation: This independence assumption applies to variance of each origin period.
End of explanation
"""


(mack.total_process_risk_**2 + mack.total_parameter_risk_**2).to_frame().values[0, -1] == \
(mack.total_mack_std_err_**2).values[0,0]

"""
Explanation: Independence is also assumed to apply to the overall standard error of reserves, total_mack_std_err_.
End of explanation
"""


mack.mack_std_err_[mack.mack_std_err_.development==mack.mack_std_err_.development.max()]

"""
Explanation: This over-reliance on independence is one of the weaknesses of the MackChainladder method. Nevertheless, if the data align with this assumption, then total_mack_std_err_ is a reasonable esimator of reserve variability.
Mack Reserve Variability
The mack_std_err_ at ultimate is the reserve variability for each origin period.
End of explanation
"""


mack.summary_

plot_data = mack.summary_.to_frame()
g = plot_data[['Latest', 'IBNR']] \
    .plot(kind='bar', stacked=True,
          yerr=pd.DataFrame({'latest': plot_data['Mack Std Err']*0,
                             'IBNR': plot_data['Mack Std Err']}),
          ylim=(0, None), title='Mack Chainladder Ultimate')
g.set_xlabel('Accident Year')
g.set_ylabel('Loss');

dist = pd.Series(np.random.normal(mack.ibnr_.sum(),
                           mack.total_mack_std_err_.values[0, 0], size=10000))
dist.plot(
    kind='hist', bins=50,
    title="Normally distributed IBNR estimate with a mean of " + '{:,}'.format(round(mack.ibnr_.sum(),0))[:-2]);

"""
Explanation: These are probably easier to see in the summary_ of the MackChainladder model.
End of explanation
"""


samples = cl.BootstrapODPSample(n_sims=10000).fit(tri['CumPaidLoss']).resampled_triangles_

"""
Explanation: ODP Bootstrap Model
The MackChainladder focused on a regression framework for determining the variability of reserve estimates.  An alternative approach is to use statistical bootstrapping or sampling from a triangle with replacement to simulate new triangles.
Bootstrapping imposes less model constraints than the MackChainladder which allows for greater applicability in different scenarios.  Sampling new triangles can be accomplished through the BootstrapODPSample estimator.  This estimator will take a single triangle and simulate new ones from it.
Notice how easy it is to simulate 10,000 new triangles from an existing triangle by accessing the resampled_triangles_ attribute.
End of explanation
"""


samples = cl.BootstrapODPSample(n_sims=10000).fit_transform(tri['CumPaidLoss'])

"""
Explanation: Alternatively, we could use BootstrapODPSample to transform our triangle into a resampled set.
End of explanation
"""


difference = round(1 - cl.Chainladder().fit(samples).ibnr_.sum('origin').mean() / \
                       cl.Chainladder().fit(tri['CumPaidLoss']).ibnr_.sum())
print("Percentage difference in estimate using original triangle and BootstrapODPSample is " +str(difference))

"""
Explanation: The notion of the ODP Bootstrap is that as our simulations approach infinity, we should expect our mean simulation to converge on the basic Chainladder estimate of of reserves.
Let's apply the basic chainladder to our original triangle and also to our simulated triangles to see whether this holds true.
End of explanation
"""


samples

pipe = cl.Pipeline([
    ('dev', cl.Development(average='simple')),
    ('tail', cl.TailConstant(1.05))])
pipe.fit(samples)

"""
Explanation: Using deterministic methods with Bootstrap samples
Our samples is just another triangle object with all the functionality of a regular triangle.  This means we can apply any functionality we want to our samples including any deterministic methods we learned about previously.
End of explanation
"""


pipe.named_steps.dev.cdf_

"""
Explanation: Now instead of a single cdf_ vector, we have 10,000.
End of explanation
"""


orig_dev = cl.Development(average='simple').fit(tri['CumPaidLoss'])
resampled_ldf = pipe.named_steps.dev.ldf_
print("12-24 LDF of original Triangle: " + str(round(orig_dev.ldf_.values[0,0,0,0],4)))
pd.Series(resampled_ldf.values[:, 0, 0, 0]).plot(
    kind='hist', bins=100,
    title='Age 12-14 LDF distribution using Bootstrap');

"""
Explanation: This allows us to look at the varibility of any fitted property used in our prior tutorials.
End of explanation
"""


mack_vs_bs = resampled_ldf.std('index').to_frame().append(
    orig_dev.std_err_.to_frame()).T
mack_vs_bs.columns = ['Mack', 'Bootstrap']
mack_vs_bs.plot(kind='bar', title='Mack Regression Framework LDF Std Err\nvs\nBootstrap Simulated LDF Std Err');

"""
Explanation: Comparison between Bootstrap and Mack
We should even be able to approximate some of the Mack parameters calculated using the regression framework.
End of explanation
"""


ibnr = cl.Chainladder().fit(samples).ibnr_.sum('origin')
ibnr_99 = ibnr.quantile(q=0.99)
print("99%-ile of reserve estimate is " +'{:0,}'.format(round(ibnr_99,0)))

"""
Explanation: While the MackChainladder produces statistics about the mean and variance of reserve estimates, those have to be fit to a distribution using MLE, MoM, etc to see the range of outcomes of reserves.  With BootstrapODPSample based fits, we can use the empirical distribution directly if we choose to.
End of explanation
"""


ax = ibnr.plot(kind='hist', bins=50, alpha=0.7, color='green').plot()
dist.plot(kind='hist', bins=50, alpha=0.4, color='blue', title='Mack vs Bootstrap Variability');

"""
Explanation: Let's see how the MackChainladder reserve distribution compares to the BootstrapODPSample reserve distribution.
End of explanation
"""


tri['EarnedPremNet'].latest_diagonal

"""
Explanation: Expected loss methods with Bootstrap
So far, we've only applied the multiplicative methods (i.e. basic chainladder) in a stochastic context.  It is possible to use an expected loss method like the BornhuetterFerguson. 
To do this, we will need an exposure vector.
End of explanation
"""


import chainladder as cl
import numpy as np
clrd = cl.load_sample('clrd')
np.prod(clrd, axis=3)

import numpy as np

np.tan(cl.load_sample('raa')).T.plot()

"""
Explanation: Passing an apriori_sigma to the BornhuetterFerguson estimator tells it to consider the apriori selection itself as a random variable.  Fitting a stochastic BornhuetterFerguson looks very much like the determinsitic version.
End of explanation
"""


# Grab completed triangle replacing simulated known data with actual known data
full_triangle = bf.full_triangle_ - bf.X_ + tri['CumPaidLoss']
# Limiting to the current year for plotting
current_year = full_triangle[full_triangle.origin==full_triangle.origin.max()].to_frame().T

"""
Explanation: We can use our knowledge of Triangle manipulation to grab most things we would want out of our model.
End of explanation
"""


# Plot the data
current_year.iloc[:, :200].reset_index(drop=True).plot(
    color='green', legend=False, alpha=0.1,
    title='Current Accident Year Expected Development Distribution', grid=True);

"""
Explanation: As expected, plotting the expected development of our full triangle over time from the Bootstrap BornhuetterFerguson model fans out to greater uncertainty the farther we get from our valuation date.
End of explanation
"""


