printing the hessian from Maxlik (Gauss 10)

Hi!

I'm getting "Hessian Calculation failed".

I have read the forum but I still cannot print it.

This is my code:

(Thanks in advance)

library maxlik;
maxset;

/* clear all globals */
clearg b, L;

/****************************************************************************************************************
modification area for different runs
********************************************************************************************************/

cls;

/* provide complete path to GAUSS dataset */
dataset = xlsreadm("C:\gauss10\\TUData3","a3",1,"");
nobs = rows(dataset);

/* pulling out the appropriate columns from the dataset */

//Aca va el numero de las columnas donde estan los TIEMPOS de las actividades que pertenecen a cada conjunto
let tAfindx = {20 26 27}; // activity free // leis, other, sleep
let tArindx = {18 19 21 23 24}; // activity restricted // Travel, Pcare, Educ, help, admin
let tADomownindx = {22 25}; // activity own domestic // hhchores, child
let tADomhiredindx = {40 41}; // activity hired domestic // random hhchores, random child

//Aca va el numero de las columnas donde estan los GASTOS de las actividades que pertenecen a cada conjunto
let eGfindx = {33 39 39}; // gasto free - leisure // leis, other, sleep
let eGrindx = {28 30 31 39 45}; // gasto restricted // Travel, Pcare, Educ, help, admin
let eDomownindx = {29 32}; // gasto own domestic // hhchores, child
let eDomhiredindx = {42 43}; // gasto hired domestic // 0.3 tot-all e, 0.15 tot-all e

//columna donde esta Tw
let tworkindx = {17};

//columna donde esta la tasa salarial, ingreso fijo y gastos fijos
let windx = {16};
let FIindx = {36};
let cfindx = {34};

L= rows(tAfindx');
P= rows(tADomownindx');
Avail = dataset[.,38]~dataset[.,38]~dataset[.,38]; // availability of alternatives
Availdom = dataset[.,38]~(dataset[.,3].>0); // availability of alternatives domestic
Tc = dataset[.,tArindx]; // committed time
Ec = sumr(dataset[., eGrindx~cfindx]) - dataset[.,FIindx]; // committed expenses

//Aca puse un tiempo minimo de recreacion (era para no dividir por cero)
treshold = 0.3; // in hours
TAf = treshold*ones(nobs,L);
TAf = substute(dataset[.,tAfindx], dataset[.,tAfindx].0))*ones(nobs,1);
i=i+1;
endo;
sigmasdom = substute(sigmasdom, sigmasdom.==0, sigmaAvdom);
sigmasdom = sigmasdom .* Availdom;

// Starting point
b = {
1.239
1.9
1.5242
1.70
1.3
1.3
1.2994
};

b = b';

/* enter title for the model */
_title = "Time use model - LISS database ";
_max_ParNames = "T_work"|"T_leis"|"T_other"|"T_sleep"|"T_chores"|"T_child"|"Psi";

/* *******************************************************************************************
MAIN PROGRAM AREA
********************************************************************************************/

//_max_GradProc = &lgd;
_max_Options = { bfgs stepbt };
_max_CovPar = 1;
_max_active = ones(L + P + 1,1)|1;
//_max_active = ones(L + 1,1)|1|0|1|1|1|0|0;

{x,f,g,cov,retcode} = maxprt(maxlik(dataset,0,&lpr,b));

clearg _max_gradproc;
b1 = x[1]|exp(x[2:L+1])|exp(x[L+P:L+P+1])|(exp(x[L+P+2]));
_max_Maxiters = 0;
{x1,f,g,cov,retcode } = maxprt(maxlik(dataset,0,&lpr1,b1));

//*************** Value of time
vtime = VTR(x1[L+P+2], dataset);
vwork = VW(vtime, dataset);

/*print("Average (Std. dev.) of Value of time as a resource [euros/hr]:", meanc(vtime)~stdc(vtime));
print("Average (Std. dev.) of Value of work [euros/hr]:", meanc(vwork)~stdc(vwork));
print("Average (Std. dev.) of Wage rate [euros/hr]:", meanc(dataset[.,windx])~stdc(dataset[.,windx]));
print("Average (Std. dev.) of Leisure time - wage rate ratio:", meanc(vtime ./ dataset[.,windx])~stdc(vtime ./ dataset[.,windx]));
print("Average (Std. dev.) of Work time - wage rate ratio:", meanc(vwork ./ dataset[.,windx])~stdc(vwork ./ dataset[.,windx]));*/

print(meanc(vtime)~stdc(vtime));
print(meanc(vwork)~stdc(vwork));
print(meanc(dataset[.,windx])~stdc(dataset[.,windx]));
print(meanc(vtime ./ dataset[.,windx])~stdc(vtime ./ dataset[.,windx]));
print(meanc(vwork ./ dataset[.,windx])~stdc(vwork ./ dataset[.,windx]));

/************************ likelihood function**********************************/
proc lpr(x,dta);
local Jac, phi, theta, Fvector, i, j, k, count, h, w, Tk, Tdomk, nk, ndomk, at, ae, sigmaTk, sigmadomTk, thetaTk, thetadomTk, timedomTk, timeTk, C, D, z, lnz;

phi = exp(x[L+P+2]);
theta = x[1]|exp(x[2:L+1])|exp(x[L+P:L+P+1])|(exp(x[L+P+2]));

// Terms to use later
ae = 1./(dta[.,windx] .* dta[.,tworkindx] - Ec - sumr(dta[.,eGfindx]) - sumr(dta[.,edomownindx]) - sumr(dta[.,edomhiredindx]));
at = sumr(dta[.,tAfindx]);
at = 1./(substute(at, at.==0, treshold));

// Jacobian
C = phi * (ae .^ 2);
D = at .^ 2;

z = zeros(nobs,1);
i = 1;
do until i == nobs+1;
Tk = (dta[i,tAfindx] .!= 0);
nk = sumr(Tk);
thetaTk = selif(theta[2:L+1], Tk'.==1);
timeTk = selif(dta[i,tAfindx]', Tk'.==1);

Tdomk = (dta[i,tADomownindx] .!= 0);
ndomk = sumr(Tdomk);
thetadomTk = selif(theta[L+P:L+P+1], Tdomk'.==1);
timedomTk = selif(dta[i,tADomownindx]', Tdomk'.==1);
sigmaTk = selif(sigmas[i,.]', Tdomk'.==1);
sigmadomTk = selif(sigmasdom[i,.]', Tdomk'.==1);

Jac = ones(nk + ndomk + 1, nk + ndomk + 1);
Jac[1,1] = -theta[1] * (dta[i,tworkindx]^(-2)) - (dta[i,windx]^2) * C[i] - D[i];

h = 2;
do until h == nk + 2;
Jac[1,h] = 0;
Jac[h,1] = -D[i];

k = 2;
do until k == nk + 2;
Jac[k,h] = -(h == k) * thetaTk[k-1] * (dta[i,tAfindx[k-1]]^(-2));
k = k + 1;
endo;
k = nk + 2;
do until k == nk + ndomk + 2;
Jac[k,h] = 0;
Jac[h,k] = -D[i];
k = k + 1;
endo;
h = h + 1;
endo;
h = nk + 2;
do until h == nk + ndomk + 2;
Jac[1,h] = dta[i,windx] * sigmadomTk[h-(nk+1)] * C[i] - D[i];
Jac[h,1] = -dta[i,windx] * sigmaTk[h-(nk+1)] * C[i] - D[i];

k = nk + 2;
do until k == nk + ndomk + 2;
Jac[k,h] = (-(h == k) * thetadomTk[k-(nk+1)] * (dta[i,tAdomownindx[k-(nk+1)]]^(-2)))+ (sigmaTk[h-(nk+1)]) * (sigmadomTk[h-(nk+1)]) * C[i] - D[i];
k = k + 1;
endo;
h = h + 1;
endo;

// Density function
Fvector = zeros(L+P+1, 1);
Fvector[1] = theta[1] / dta[i,tworkindx] + phi * dta[i,windx] * ae[i] - at[i];

j = 2;
do until j == L + 2;
if (Tk[j-1] == 1);
Fvector[j] = theta[j] / dta[i,tAfindx[j-1]] - at[i];
endif;
j = j + 1;
endo;

j = L + 2;
do until j == L + P + 2;
if (Tdomk[j-(L+1)] == 1);
Fvector[j] = theta[j] / dta[i,tADomownindx[j-(L+1)]] - at[i] - phi * sigmas[i,j-(L+1)]*ae[i];
endif;
j = j + 1;
endo;

z[i] = abs(det(J)) * pdfmvn(Fvector, 0, eye(L+P+1));
i = i + 1;
endo;

w = zeros(rows(dta),1);

if z > w;
lnz = ln(z);
else;
lnz = ln(z-((z. w;
lnz = ln(z);
else;
lnz = ln(z-((z.<=w).*(z-0.0001))); endif; retp(lnz); endp; /********************** end of likelihood function*****************************/ /************************ PDF Normal **********************************/ proc pdfmvn(x,u,s); local d,p1,p2; d = rows(x); p1 = exp(-0.5*((x-u)'*inv(s)*(x-u))); p2 = ((2*pi)^(d/2))*sqrt(det(s)); retp(p1/p2); endp; /************************ Value of time as resource **********************************/ proc VTR(psi,dta); local ae1, at, vtr; ae1 = (dta[.,windx] .* dta[.,tworkindx] - Ec - sumr(dta[.,eGfindx])- sumr(dta[.,edomownindx])- sumr(dta[.,edomhiredindx])); at = sumr(dta[.,tAfindx]); at = 1./(substute(at, at.==0, treshold)); retp(psi .* at .* ae1); endp; proc VW(vtr, dta); retp(vtr - dta[.,windx]); endp;

4 Answers

0

In Maxlik, the Hessian used for computing the covariance matrix of the parameters is stored in a global, _max_FinalHess.  When the Hessian fails to invert, put the statement

print _max_FinalHess;

after the call to Maxlik.

The failure to invert generally means the Hessian contains linear dependencies.  One reason might be that some parameters are not identified, i.e., there isn't enough information in the data to produce an estimate.  It also could be a scaling issue.  Be sure that your data is scaled properly.  Mixing very large numbers with small numbers can be computationally catastrophic.

A method to find the linear dependencies is described in the final section of a white paper, Optimization with the Quasi-Newton Method.

First compute

{ r,e } = qre(_max_finalHess);

Zeros on the diagonal of r indicate linear dependencies.  The e vector indicates the rows of _max_finalHess associated with the rows of r.  Suppose the last two diagonal elements of r are essentially zero, then the last two elements of e indicate the linearly dependent rows of _max_finalHess.

The description of the linear dependencies can be found by computing

b = utrisol(r[1:rows(r)-2,rows(r)-1:rows(r)],r[1:rows(r)-2,1:rows(r)-2))

b will have two columns describing the exact linear dependency of the rows of _max_finalHess as indicated by the last two elements of e, to the remaining rows of _max_finalHess.

link

RonSchoenberg

210

0

Thanks Ron!!!!

I got:

    0.0049   -4.3987e-6    -3.25011e-7    -1.2195e-5   0.0199   0  0 
-4.3987e-6        0.032     -6.5749e-7    -2.8294e-7  -0.0114   0  0
-3.2501e-7   -6.5749e-7         0.0067     -2.017e-7  -0.0017   0  0 
 -1.219e-5    -2.829e-7      -2.017e-7        0.0003 -1.27e-7   0  0 
    0.0199      -0.0114        -0.0017     -1.272e-7   0.4474   0  0
         0            0              0             0        0   0  0
         0            0              0             0        0   0  0

Everything = 0 (almost).

What can I do?

Thankssss again!
And have a great 2015.

link

jorgerosales

0

0

The Hessian is a 7X7 matrix.

🙂

link

jorgerosales

0

0

It appears that you have two linear dependencies.  Those last two rows are clearly zero.

The Hessian is very nearly entirely zero though, suggesting that your likelihood function is basically flat.  You either have a serious data problem, or a serious scaling problem.

link

RonSchoenberg

210