ARIMA estimation with ML

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 estimation of sigma and mu in the code at the end?

Note: This is already written in the comment of the thread : Mathworks, however as it is a new topic I though it would be better to open up another thread.

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 :)