%% example 4.3.1
iRate = .04
v = 1/(1+iRate)

iltable % loads data from table

A(40) = 0; % initial value, A^1_40:0

for k = 9:-1:0
    A(30+k) = v*q(30+k) + v*(1-q(30+k))*A(30+k+1);
end
disp(A(30:40))
disp(['val for insurance ' num2str(A(30))])

A2(40) = 0; % initial value, 2A^1_40:0

for k = 9:-1:0
    A2(30+k) = v^2*q(30+k) + v^2*(1-q(30+k))*A2(30+k+1);
end

disp(['second moment ' num2str(A2(30))])
VarianceA = A2(30) - (A(30))^2

%% continuation -- suppose 
% 100 policies issued
% each pays $1000
%
% S = portfolio
%
PMT = 1000
nPolicies = 100
ExpectedS = PMT*nPolicies*A(30)
VarianceS = nPolicies*PMT^2*VarianceA
% find inverse of c.d.f. of Normal(0,1)
NI = norminv(.95, 0, 1)
h = NI*sqrt(VarianceS) + ExpectedS
disp(['fund should start with ' num2str(h)])
disp(['compare with expected value ' num2str(ExpectedS)])