3dplot.R

library(plotly)

## Loading required package: ggplot2
## 
## Attaching package: 'plotly'
## 
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## 
## The following object is masked from 'package:graphics':
## 
##     layout

## Loading required package: ggplot2
## 
## Attaching package: 'plotly'
## 
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## 
## The following object is masked from 'package:graphics':
## 
##     layout
# Life Table

data = read.csv("C:/Users/Nathan/Downloads/ca80m.csv") 

# Parameters
delta  = 0.05
delta0 = 0.03
alpha  = 0.1
sigma  = 0.02


# Mth moment of conditional present value function/discount function
# When force of interest is modeled by Ornstein-Uhlenbeck process
#
# @param k curtate future lifetime
# @param m moment

EVDiscount <- function(k, m) {
    
    EV = ((delta0-delta)/alpha)*(1-exp(-alpha*(k+1))) + (k+1)*delta
    V  = ((sigma/alpha)^2)*(k+1) + (sigma^2)/(2*alpha^3)*(-3 + 4*exp(-alpha*(k+1)) - exp(-2*alpha*(k+1)))
    
    exp(-m*EV + (m^2)*0.5*V)
}

# This is a tool to calculate mth moment of present value of the benefit
# 
# @param x issue age
# @param n term
# @param m is moment 
#
# If n = omega-x, this is a whole life insurance

f <- function(x, n, m=1) {
    
    # Generate a table of k|qx, k = 0, ..., 10
    kpx = rep(1, n+1)
    for (k in 1:n) {
        kpx[k+1] = kpx[k] * data[x+k,3]
    }
    
    qx = numeric(n+1)
    for (k in 0:(n-1)) {
        qx[k+1] = data[x+k+1,2]
    }
    
    kqx = numeric(n+1)
    for (k in 0:(n-1)) {
        kqx[k+1] = kpx[k+1]* qx[k+1]
    }
    
    table = data.frame(k=0:n, kpx=kpx, qx=qx, kqx=kqx)
    
    k = 0:(n-1)
    vk = sapply(k, function(x) EVDiscount(x,m))
    
    A <- function(x) table[x+1,4]*vk[x+1]
    sum(sapply(k, A))
}

g <- function(x,n) {
    sqrt(f(x,n,2) - f(x,n,1)^2)
}

x  = 20:70
n  = 1:50
nx = length(x)
ny = length(n)
z  = matrix(0, nx, ny)

for(i in 1:nx) {
    for(j in 1:ny) {
        z[i,j] = g(x[i],n[j])
    }
}

plot_ly(x=x,y=n,z=z,type="surface",text=labels)