A Maximum Likelihood estimation with fminunc
Show older comments
Hello,
I would like to estimate the parameters k1 and k2 using the maximum likelihood method. To check the result, I have previously made an estimate using the OLS method (OLSk1 and OLSk2).
The problem now is that although k2 is identical in both models, k1 is not. Since OLSk1 is close to values from papers, the estimation of the maximum likelihood method must be wrong, but I don't get an error. Is perhaps the 'quasi-newton' algorithm the wrong choice for a log-function?
Thanks in advance for your help.
Kind regards :)
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Here the OLS estimation, further down the labeled maximum likelihood.
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
[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)';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%here does the Maximum Likelihood estimation starts
results = zeros(n, 2);
for i = CellyearStart:CellyearEnd
iter = 10^10;
dura = 0;
x0=[-2, 0.1];
fun = @(k) sum_l(k, (i-1)*globalageEnd+1, (i)*globalageEnd, EInitial_xt, D_xt, E_xt);
options=optimoptions(@fminunc,'Algorithm','quasi-newton');
while iter>0.00000001 && dura<4000000
dura = dura + 1;
x1 = fminunc(fun, x0, options);
iter = sum(abs(x1-x0));
x0 = x1;
end
results(i-CellyearStart+1, :) = x1;
end
% plot(time, results(:, 1), 'r')
% hold on
% plot(time, kt1, 'b')
function ml = l(k, x, D, E)
k1 = k(1);
k2 = k(2);
ml = (-1)*(D*log(E*log(1+exp(k1+(x-52.5)*k2)))-E*log(1+exp(k1+(x-52.5)*k2))-log(sqrt(2*pi*D))-D*log(D/exp(1)));
end
function sum_ml = sum_l(k, startIndex, endIndex, EInitial_xt, D_xt, E_xt)
sum_ml = 0 ;
for j = startIndex:endIndex
sum_ml = sum_ml + l(k, EInitial_xt(j, 2), D_xt(j, 1), E_xt(j, 1));
end
end
1 Comment
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 :)
Accepted Answer
By the way: Why do you iterate the solutions obtained from fminunc ?
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Here the OLS estimation, further down the labeled maximum likelihood.
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
[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)';
k1 = -5.3:0.1:-3.3;
k2 = 0.04:0.01:0.12;
iyear = CellyearStart;
for i = 1:numel(k1)
for j = 1:numel(k2)
F(i,j) = sum_l([k1(i),k2(j)], (iyear-1)*globalageEnd+1, (iyear)*globalageEnd, EInitial_xt, D_xt, E_xt);
end
end
surf(k2,k1,F)
function ml = l(k, x, D, E)
k1 = k(1);
k2 = k(2);
ml = (-1)*(D*log(E*log(1+exp(k1+(x-52.5)*k2)))-E*log(1+exp(k1+(x-52.5)*k2))-log(sqrt(2*pi*D))-D*log(D/exp(1)));
end
function sum_ml = sum_l(k, startIndex, endIndex, EInitial_xt, D_xt, E_xt)
sum_ml = 0 ;
for j = startIndex:endIndex
sum_ml = sum_ml + l(k, EInitial_xt(j, 2), D_xt(j, 1), E_xt(j, 1));
end
end
More Answers (1)
Is perhaps the 'quasi-newton' algorithm the wrong choice for a log-function?
More likely, the solution is correct, but your expectations are not.
Since you only have a 2 parameter model, you can easily plot the loglikelihood as a surface to see if the correct minima has been found.
Categories
Find more on Solver Outputs and Iterative Display in Help Center and File Exchange
