Hello,
context of the Cairns Blake Dowd model (2006):
I want to estimate the paramteres of the bivariate random walk with drift, i.e., mu and V. Where V is the variance-covariance-matrix and mu is the drift component of the random walk.
We need V in order to use the Cholesky decomposition to calculate C.
I estimated the mu and sigma due to the OLS and the covariance function in Matlab, both values are the quite same. However, they differ from the values in the paper, which was also the case for the k1 estimation from our maximum likelihood perviously in the code (I think the paper might be wrong here)
And if we compare the rest of the values, for the mu and sigma from our maximum likelihood at the bottom of the code, the values are quite the same, despite the value for the k2 sigma. Nevertheless, the comparision between mu and sigma from the OLS estimation and the maximum likelihood estimation reveals that these values are substantially different.
Is something wrong with the code?
EInitial_xt = xlsread('Exposures.xls');
E_xt = EInitial_xt(: ,4); %populatiion by age and year
QInitial_xt = xlsread('qMale.xlsx');
Q_xt = QInitial_xt(: ,4);
DInitial_xt = xlsread('Deaths.xls');
D_xt = DInitial_xt(: ,4); %Death year, age
globalyearStart = 1841;
globalyearEnd = 2021;
globalageStart = 0;
globalageEnd = 105;
yearStart = 1960;
yearEnd = 2006;
ageStart = 60;
ageEnd = 89;
m0 = globalageEnd - globalageStart + 1; % quantity of total ages
n0 = globalyearEnd - globalyearStart + 1; % quantity of total years
m = ageEnd - ageStart+1; % quantity of age considered
n = yearEnd - yearStart+1; % quantity of years considered
x = zeros(m, 1);
for i= 1:m
x(i, 1) = i-1+ageStart;
end
globalorganizedtableE = zeros(m0+1, n0+1); % all data E_xt
globalorganizedtableD = zeros(m0+1, n0+1);
globalorganizedtableQ = zeros(m0+1, n0+1);
for k = 1:m0*n0
age = mod(k-1, m0) + 1; % cyclical arithmetic
yearIndex = ceil(k/ m0); % Ceiling in order to skip through the columns after the span of ages m
globalorganizedtableE(age+1, yearIndex+1) = EInitial_xt(k, 4);
globalorganizedtableD(age+1, yearIndex+1) = DInitial_xt(k, 4);
globalorganizedtableQ(age+1, yearIndex+1) = QInitial_xt(k, 4);
end
for i = globalageStart+1:globalageEnd+1
globalorganizedtableE(i+1, 1)= i-1;
globalorganizedtableD(i+1, 1)= i-1;
globalorganizedtableQ(i+1, 1)= i-1;
end
for j = 1:n0
globalorganizedtableE(1, j+1) = EInitial_xt(m0*j, 1);
globalorganizedtableD(1, j+1) = EInitial_xt(m0*j, 1);
globalorganizedtableQ(1, j+1) = EInitial_xt(m0*j, 1);
end
% structured tables of the total values only from the considered observations.
organizedtableE = globalorganizedtableE((ageStart+2):(ageEnd+2), (yearStart-globalyearStart+2):(yearEnd-globalyearStart+2));
organizedtableD = globalorganizedtableD((ageStart+2):(ageEnd+2), (yearStart-globalyearStart+2):(yearEnd-globalyearStart+2));
organizedtableQ = globalorganizedtableQ((ageStart+2):(ageEnd+2), (yearStart-globalyearStart+2):(yearEnd-globalyearStart+2));
CellyearStart= yearStart -1841+1; % wanted starting year - start year of data +1
CellyearEnd= yearEnd -1841+1; % wanted ending year - start year of data +1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% OLS estimation
[a, b] = size(organizedtableQ);
logitq = log(organizedtableQ./(ones(a,b)-organizedtableQ));
sumlq = sum(logitq);
meanx= mean(x);
xf=x-meanx;
xf2 = xf.^2;
sumxf = sum(xf);
sumxf2 = sum(xf2);
num2 = (length(x).*xf')*logitq;
den2 = length(x)*(sumxf2);
OLSk2 = num2./den2;
OLSk1 = (sumlq)./length(x);
time =(yearStart:yearEnd)';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Maximum Likelihood estimtaion for k1 and k2
results = zeros(n, 2);
for i = CellyearStart:CellyearEnd
x0=[0, 0.2];
fun = @(k) sum_loglikelihood(k, (i-1)*globalageEnd+1, (i)*globalageEnd, EInitial_xt, D_xt, E_xt, meanx);
options=optimoptions(@fminunc,'Algorithm','quasi-newton');
x1 = fminunc(fun, x0, options);
results(i-CellyearStart+1, :) = x1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Arima (0, 1, 0), bivariat random walk with drift
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% OLS estimation
kt = [OLSk1, OLSk2];
a = (kt(:,n) - kt(:,1))/(yearEnd - yearStart + 1);
ktt1 = kt(:,1:(end-1));
ktt2 = kt(:,2:end);
Dkt = ktt2-ktt1-a*ones(1,yearEnd - yearStart);
sigma = sum(Dkt.^2,2)/(yearEnd - yearStart + 1);
% Cholesky decomposition of covariance matrix
Corr = corrcoef(OLSk1, OLSk2);
C = chol(Corr);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Maximum likelihood estimation for Arima (0, 1, 0) model.
resultsARIMA = zeros(2, 2);
for i = 1:2
x0=[-0.0159028602730627, 0.000982643838986423];
funARIMA = @(SigmaMu) sum_loglikelihood_ARIMA (SigmaMu, n, i, Delta);
options = optimoptions(@fminunc,'Algorithm','quasi-newton');
x1 = fminunc(funARIMA, x0, options);
resultsARIMA(i, :) = x1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% functions
function ml = loglikelihood(k, x, D, E, meanx)
k1 = k(1);
k2 = k(2);
ml = (-1)*(D*log(E*log(1+exp(k1+(x-meanx)*k2)))-E*log(1+exp(k1+(x-meanx)*k2))-log(sqrt(2*pi*D))-D*log(D/exp(1)));
end
function sum_ml = sum_loglikelihood(k, startIndex, endIndex, EInitial_xt, D_xt, E_xt, meanx)
sum_ml = 0 ;
for j = startIndex:endIndex
sum_ml = sum_ml + loglikelihood(k, EInitial_xt(j, 2), D_xt(j, 1), E_xt(j, 1), meanx);
end
end
function ml_ARIMA = loglikelihood_ARIMA(SigmaMu ,Delta)
mu = SigmaMu(1);
sigma = SigmaMu(2);
ml_ARIMA = log(sigma*sqrt(2*pi)) + 0.5*((Delta-mu)/sigma)^2;
end
function sum_ml_ARIMA = sum_loglikelihood_ARIMA(SigmaMu, DurationYear, RowIndex, Delta)
sum_ml_ARIMA = 0;
for i = 1:DurationYear-1
sum_ml_ARIMA = sum_ml_ARIMA + loglikelihood_ARIMA(SigmaMu, Delta(RowIndex, i));
end
end
Kind regards :)
