Error G0058 : Index out of range in likelihood procedure

Dear all,
I am getting:

C:\gauss10\LISSModel_all_iidxxx(289) : error G0058 : Index out of range
Currently active call: lpr [289] C:\gauss10\LISSModel_all_iidxxx
Stack trace:
lpr called from C:\gauss10\maxutil.src, line 2910
_max_rdd called from C:\gauss10\maxutil.src, line 405
_max called from C:\gauss10\maxlik.src, line 515
maxlik called from C:\gauss10\LISSModel_all_iidxxx, line 108

When running this code:

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 = {
15.239
0.9
2.5242
4.70
1.3
1.3
-2.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-1] * C[i] - D[i];
Jac[h,1] = -dta[i,windx] * sigmaTk[h-1] * C[i] - D[i];

k = nk + 2;
do until k == nk + ndomk + 2;
Jac[k,h] = (-(h == k) * thetadomTk[k-1] * (dta[i,tADomownindx[k-1]]^(-2))) + (sigmaTk[h-1]) * (sigmadomTk[h-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 = 5;
do until j == L+P+2;
if (Tdomk[j-1] == 1);
Fvector[j] = theta[j] / dta[i,tADomownindx[j-1]] - at[i] - phi * sigmas[i,j-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;

6 Answers

0

This error message is telling us that the error is on line 289 of your file LISSModel_all_iidxxx, inside of the procedure lpr. The error "Index out of range", means most likely that your program is trying to reference a row or column that does not exist. For example, if you had a matrix with 3 rows and 2 columns and tried to index the element from the 4th row and 2nd column, you would get the error "Index out of range":

//create a 3x2 matrix
x = { 2 4,
      5 1,
      8 3 };
//the line below will return index out of range,
//because there is no 4th row in 'x'
y = x[4, 2];

You can either use the debugger or place a print statement before the line that is causing the problem to see which index is causing the problem.

link

aptech

1,773

0

Thanks!!!!
I will work on that.

link

jorgerosales

0

0

This is the line:

Jac[1,h] = dta[i,windx] * sigmadomTk[h-1] * C[i] - D[i];

But I'm not requesting something out of range.

🙁

link

jorgerosales

0

0

Try adding these lines and then running your program to see what you find:

print "************************";
print "h = " h;
print "cols(Jac) = " cols(Jac);
print "-------------------------";
print "rows(dta) = " rows(dta);
print "rows(C) = " rows(C);
print "rows(D) = " rows(D);
print "-------------------------";
print "cols(dta) = " cols(dta);
print "windx = " windx;
print "i = " i;
print "-------------------------";
print "rows(sigmadomTk) = " rows(sigmadomTk);
print "h - 1 = " (h - 1);

Jac[1,h] = dta[i,windx] * sigmadomTk[h-1] * C[i] - D[i];

I think you will find the indexing problem in the final iteration of these print statements. If not, post the last output that you see from the program. (NOTE: I made a guess that C and D were column vectors. If I am wrong, you will need to change "rows(C)" to "cols(C)")

link

aptech

1,773

0

This is the last output

What should I look for?

h = 2.0000000
cols(Jac) = 3.0000000
-------------------------
rows(dta) = 1193.0000
rows(C) = 1193.0000
rows(D) = 1193.0000
-------------------------
cols(dta) = 47.000000
windx = 16.000000
i = 133.00000
-------------------------
rows(sigmadomTk) = 1.0000000
h - 1 = 1.0000000

link

jorgerosales

0

0

I fixed it, but I cannot print the Hessian.

Can you help me please?

link

jorgerosales

0