ACTL3143 & ACTL5111 Deep Learning for Actuaries
Customer 1 = (25, 3, 🚙)
Customer 2 = (40, 5, 🚐)
Customer 3 = (19, 1, 🏎️)
Customer 4 = (60, 10, 🚘)
All customers
# Plot all on the same graph
plt.figure()
for i, (age, car_age, car_type, cost) in enumerate(data):
mean = cost
shape = 2 # Example shape parameter for the gamma distribution
scale = mean / shape # Scale parameter to achieve the correct mean
plt.plot(x, stats.gamma.pdf(x, a=shape, scale=scale), label=f"Customer {i+1}", color=colors[i])
plt.title(f"Predicted Claim Size Distribution")
plt.xlabel("Claim size")
plt.ylabel("Density")
plt.legend()Traditional:
| 👤 | 🚗 | 🏎️ |
|---|---|---|
| 35 | 12 | Sport |
Model→
1582.48
Distributional:
| 👤 | 🚗 | 🏎️ |
|---|---|---|
| 35 | 12 | Sport |
Model→
ages = np.linspace(18, 80, 200) # Driver ages from 18 to 80
claim_sizes = 1000 * (1 - np.exp(-((ages - 35) / 20)**2))
# Create the updated plot
# plt.figure(figsize=(8, 6))
plt.plot(ages, claim_sizes, linewidth=2)
plt.title("Effect of Driver Age on Claim Sizes")
plt.xlabel("Driver Age")
plt.ylabel("Claim Size Effect ($)")
# plt.grid(alpha=0.3)
# plt.axvline(35, color='gray', linestyle='--', label='New Middle Age Peak')
# plt.legend(fontsize=10)
plt.show()GLMs cannot (easily) do this \longrightarrow Use a neural network
x = np.linspace(0, 20*500, 500) # Range restricted to positive real numbers
# Adjust the components to ensure a bimodal mixture
pdf1 = stats.lognorm.pdf(x, s=0.3, scale=np.exp(1)*500) # Sharper log-normal distribution
pdf2 = stats.gamma.pdf(x, a=10, scale=0.8*500) # Narrower gamma distribution
mixture_pdf = 0.3 * pdf1 + 0.7 * pdf2 # Equal weights to emphasize bimodality
# Create the updated plot
plt.figure();
plt.plot(x, pdf1, label='Parent driving', linestyle='--', color='blue')
plt.plot(x, pdf2, label='Kid driving', linestyle='--', color='green')
plt.plot(x, mixture_pdf, label='Mixture Distribution', color='red', linewidth=2)
plt.title("Claim Size Distribution")
plt.legend()
plt.xlabel("Claim Size")
plt.ylabel("Density")
# Turn off the y axis
# plt.gca().axes.get_yaxis().set_visible(False)
None
Classifiers already give us a probability, which is a big step up compared to regression models.
However, neural networks’ “probabilities” can be overconfident.
We already saw a case of this.
See Guo et al. (2017), On Calibration of Modern Neural Networks.
# Generate sample data for linear regression
np.random.seed(0)
X_toy = np.linspace(0, 10, 10)
np.random.shuffle(X_toy)
beta_0 = 2
beta_1 = 3
y_toy = beta_0 + beta_1 * X_toy + np.random.normal(scale=2, size=X_toy.shape)
sigma_toy = 2 # Assuming a standard deviation for the normal distribution
# Fit a simple linear regression model
coefficients = np.polyfit(X_toy, y_toy, 1)
predicted_y = np.polyval(coefficients, X_toy)
# Plot the data points and the fitted line
plt.scatter(X_toy, y_toy, label='Data Points')
plt.plot(X_toy, predicted_y, color='red', label='Fitted Line')
# Draw the normal distribution bell curve sideways at each data point
for i in range(len(X_toy)):
mu = predicted_y[i]
y_values = np.linspace(mu - 4*sigma_toy, mu + 4*sigma_toy, 100)
x_values = stats.norm.pdf(y_values, mu, sigma_toy) + X_toy[i]
plt.plot(x_values, y_values, color='blue', alpha=0.5)
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.legend()
y_pred = np.polyval(coefficients, X_toy[:4])
fig, axes = plt.subplots(4, 1, figsize=(5.0, 3.0))
x_min = y_pred[:4].min() - 4*sigma_toy
x_max = y_pred[:4].max() + 4*sigma_toy
x_grid = np.linspace(x_min, x_max, 100)
# Plot each normal distribution with different means vertically
for i, ax in enumerate(axes):
mu = y_pred[i]
y_grid = stats.norm.pdf(x_grid, mu, sigma_toy)
ax.plot(x_grid, y_grid)
ax.set_ylabel(f'$f(y ; \\boldsymbol{{x}}_{{{i+1}}})$')
ax.set_xticks([y_pred[i]], labels=[r'$\mu_{' + str(i+1) + r'}$'])
ax.plot(y_toy[i], 0, 'rx', clip_on=False)
plt.tight_layout();
def lagged_timeseries(df, target, window=30):
lagged = pd.DataFrame()
for i in range(window, 0, -1):
lagged[f"T-{i}"] = df[target].shift(i)
lagged["T"] = df[target].values
return lagged
stocks = pd.read_csv("../Time-Series-And-Recurrent-Neural-Networks/aus_fin_stocks.csv")
stocks["Date"] = pd.to_datetime(stocks["Date"])
stocks = stocks.set_index("Date")
_ = stocks.pop("ASX200")
stock = stocks[["CBA"]]
stock = stock.ffill()
# Compute daily log returns
stock_log = np.log(stock / stock.shift(1)).dropna()
# Helper functions for converting log returns to prices
def log_to_price(log_returns, initial_price):
cumulative_log_returns = log_returns.cumsum()
return initial_price * np.exp(cumulative_log_returns)
def get_last_price(stock_df, cutoff_date):
last_known_date = stock_df.loc[:cutoff_date].index[-1]
return stock_df.loc[last_known_date, "CBA"]
# Create lagged features from log returns
df_lags = lagged_timeseries(stock_log, "CBA", 40)
# Split the data in time
X_train = df_lags.loc[:"2024"]
X_test = df_lags.loc["2025":]
# Remove any with NAs and split into X and y
X_train = X_train.dropna()
X_test = X_test.dropna()
y_train = X_train.pop("T")
y_test = X_test.pop("T")
lr = LinearRegression()
lr.fit(X_train, y_train);
# Plot the mean forecast
plt.figure(figsize=(8, 3))
plt.plot(stock.loc["2024-12":"2025"].index, stock.loc["2024-12":"2025"]["CBA"], label="CBA")
plt.plot(mean_forecast, label="Mean")
# Plot the quantile-based shaded area
plt.fill_between(mean_forecast.index,
lower_quantile,
upper_quantile,
color="grey", alpha=0.2)
# Plot settings
plt.axvline(pd.Timestamp("2025-01-01"), color="black", linestyle="--")
plt.legend(loc="center left", bbox_to_anchor=(1, 0.5))
plt.xlabel("Date")
plt.ylabel("Stock Price")
plt.tight_layout();
y_pred = lr.predict(X_train)
residuals = y_train - y_pred
residuals -= np.mean(residuals)
residuals /= np.std(residuals)
stats.shapiro(residuals)/home/plaub/miniconda3/envs/ai/lib/python3.11/site-packages/scipy/stats/_axis_nan_policy.py:586: UserWarning: scipy.stats.shapiro: For N > 5000, computed p-value may not be accurate. Current N is 6582.
res = hypotest_fun_out(*samples, **kwds)ShapiroResult(statistic=np.float64(0.9305463633130856), pvalue=np.float64(5.938320016175858e-48))Heteroskedasticity!
# Drop 2008-2009 and also 2020-2021
mask = ~((y_train.index.year >= 2008) & (y_train.index.year <= 2009)) & ~((y_train.index.year >= 2020) & (y_train.index.year <= 2021))
X2 = X_train.loc[mask]
y2 = y_train.loc[mask]
lr2 = LinearRegression()
lr2.fit(X2, y2)
res2 = y2 - lr2.predict(X2)
res2 = (res2 - np.mean(res2)) / np.std(res2)
plt.plot(y2.index, res2)
plt.xlabel("Date")
plt.ylabel("Standardised Residuals")
plt.tight_layout();Refit the model without these crisis years.
ShapiroResult(statistic=np.float64(0.9797545738506994), pvalue=np.float64(1.8584817750981994e-27))
X_train, X_test, y_train, y_test = train_test_split(
sev.drop("ClaimAmount", axis=1), sev["ClaimAmount"], random_state=2023)
ct = make_column_transformer((OrdinalEncoder(), ["Area", "VehGas"]),
("drop", ["VehBrand", "Region"]), remainder=StandardScaler())
X_train = ct.fit_transform(X_train)
X_test = ct.transform(X_test)
plt.hist(y_train[y_train < 5000], bins=30);
# Make some example where the distribution is multimodal because of a binary covariate which separates the means of the two distributions
np.random.seed(1)
fig, axes = plt.subplots(3, 1, figsize=(5.0, 3.0), sharex=True)
x_min = 0
x_max = y_train.max()
x_grid = np.linspace(x_min, x_max, 100)
# Simulate some data from an exponential distribution which has Pr(X < 1000) = 0.9
n = 100
p = 0.1
lambda_ = -np.log(p) / 1000
mu = 1 / lambda_
y_1 = np.random.exponential(scale=mu, size=n)
# Pick a truncated normal distribution with a mean of 1100 and std of 250 (truncated to be positive)
mu = 1100
sigma = 100
y_2 = stats.truncnorm.rvs((0 - mu) / sigma, (np.inf - mu) / sigma, loc=mu, scale=sigma, size=n)
# Combine y_1 and y_2 for the final histogram
y = np.concatenate([y_1, y_2])
# Determine common bins
bins = np.histogram_bin_edges(y, bins=30)
# Plot each normal distribution with different means vertically
for i, ax in enumerate(axes):
if i == 0:
ax.hist(y_1, bins=bins, density=True, color=colors[i+1])
ax.set_ylabel(f'$f(y | x = 1)$')
elif i == 1:
ax.hist(y_2, bins=bins, density=True, color=colors[i+1])
ax.set_ylabel(f'$f(y | x = 2)$')
else:
ax.hist(y, bins=bins, density=True)
ax.set_ylabel(f'$f(y)$')
plt.tight_layout();Combining GLM & ANN.
# Ensure reproducibility
random.seed(1)
# Make a 4x1 grid of plots
fig, axes = plt.subplots(4, 1, figsize=(5.0, 3.0), sharex=True)
# Define the x-axis
x_min = 0
x_max = 5000
x_grid = np.linspace(x_min, x_max, 100)
# Plot a few Gamma distribution pdfs with different means.
# Then plot Gamma distributions with shifted means and the same dispersion parameter.
glm_means = [1000, 3000, 2000, 4000]
cann_means = [1500, 1400, 3000, 5000]
for i, ax in enumerate(axes):
ax.plot(x_grid, stats.gamma.pdf(x_grid, a=2, scale=glm_means[i]/2), label=f'GLM')
ax.plot(x_grid, stats.gamma.pdf(x_grid, a=2, scale=cann_means[i]/2), label=f'CANN')
ax.set_ylabel(f'$f(y | x_{i+1})$')
if i == 0:
ax.legend(["GLM", "CANN"], loc="upper right", ncol=2)The CANN architecture.
An MDN that outputs the parameters for a K component mixture distribution. \boldsymbol{\theta}_k(\boldsymbol{x}; \boldsymbol{w}^*)= (\theta_{k,1}(\boldsymbol{x}; \boldsymbol{w}^*), \ldots, \theta_{k,|\boldsymbol{\theta}_k|}(\boldsymbol{x}; \boldsymbol{w}^*)) consists of the parameter estimates for the kth mixture component.
We demonstrate the structure of a Gamma MDN that outputs the parameters for a Gamma mixture with two components.
y_pred = np.polyval(coefficients, X_toy[:4])
y_pred[2] *= 1.1
sigma_preds = sigma_toy * np.array([1.0, 3.0, 0.5, 0.5])
fig, axes = plt.subplots(1, 4, figsize=(5.0, 3.0), sharey=True)
x_min = y_pred[:4].min() - 4*sigma_toy
x_max = y_pred[:4].max() + 4*sigma_toy
x_grid = np.linspace(x_min, x_max, 100)
# Plot each normal distribution with different means vertically
for i, ax in enumerate(axes):
y_grid = stats.norm.pdf(x_grid, y_pred[i], sigma_preds[i])
ax.plot(x_grid, y_grid)
ax.plot([y_toy[i], y_toy[i]], [0, stats.norm.pdf(y_toy[i], y_pred[i], sigma_preds[i])], color='red', linestyle='--')
ax.scatter([y_toy[i]], [stats.norm.pdf(y_toy[i], y_pred[i], sigma_preds[i])], color='red', zorder=10)
ax.set_title(f'$f(y ; \\boldsymbol{{x}}_{{{i+1}}})$')
ax.set_xticks([y_pred[i]], labels=[r'$\mu_{' + str(i+1) + r'}$'])
# ax.set_ylim(0, 0.25)
# Turn off the y axes
ax.yaxis.set_visible(False)
plt.tight_layout();An example of neurons dropped during training.
Combine many models to get better predictions.
