/* MixModel.cmd: This RPL model estimates nine coefficients, all fixed except a4 and a7--which are assumed random normal--and also estimates cov(a4,a7). It requires as input the data file XMAT, the file YVEC which indicates which choice was made for each choice occasion in XMAT, and the file TIMES which indicates the number of choice occasions for each individual This program was developed and made available by Kenneth Train (http://elsa.berkeley.edu/~train/software.html) and has been modified by Kathleen Greer Rossmann and Edward Morey to estimate the Mixture Model presented in "Using stated-preference questions to investigate variations in willingness to pay for preserving marble monuments: classical heterogeneity random parameters, and mixture models." To run this program, change the directory location indicated below and put XMAT, YVEC, and TIMES in the same location. */ @======================================================================@ @@@ Modify the following variables @@@ @@@ See the Train's file manual.txt for documentation. @@@ @======================================================================@ st = "C:\\gauss\\jce"; ChangeDir(st); output file=MixModel.out reset; NP = 259; NOBS = 2568; NALT = 2; load XMAT; load YVEC; load TIMES; NREP = 5000; save NREP; SEED1 = 27; save SEED1; NITER = 1000; EPS = 1.e-4; @ 1 to weight the log-likelihood function and gradient else 0 @ DOWT = 0; @ If DOWT=1, load or create an NPx1 vector of weights. @ WEIGHT = {0}; @ Number of variables in XMAT. (integer) @ NVAR =9; @ Give the number of explanatory variables that have fixed coefficients. @ @ (integer) @ NFC = 7; @ Give the number of explanatory variables that have normally @ @ distributed coefficients. (integer) @ NNC = 2; @ Give the number of explanatory variables that have log-normally @ @ distributed coefficients. (integer) @ NLC = 0; @ (NFC+NNC+NLC)x1 vector to identify the variables in XMAT that have @ @ the associated distributed coefficients. (column vector) @ IDC = { 1, 2, 3, 5, 6, 8, 9, 4, 7}; @ 1 if all people do not face all NALT alternatives, else 0. @ CENSOR = 0; @ Identify the variable in XMAT that identifies the alternatives @ @ each person faces. (integer) @ IDCENSOR = 0; @ 1 to print out the diagonal elements of the Hessian, else 0. @ PRTHESS = 1; @ 1 to rescale any variable in XMAT, else 0. @ RESCALE = 0; @ If RESCALE = 1 create a q x 2 matrix to id which variables to @ @ rescale and by what multiplicative factor, else make RESCLMAT @ @ a dummy value. @ RSCLMAT = { 0}; @ 1 if NLC=0 and you want to estimate the full var/cov matrix @ DOCH = 1; @ {NFC+(2*NNC)+(2*NLC)} x 1 vector of starting values where the first @ @ NFC+NNC+NLC are the means and the last NNC+NLC are the associated @ @ error terms. (column vector) @ @ If DOCH=1, B=NFC+NNC+sum(1+2+...+NNC) where the first NFC+NNC @ @ elements are the @ @ means and the remaining elements are the values of the corresponding @ @ Cholesky factor in row major order. @ B = {-1.2653, 0.1588, 0.0024, 1.4905, -0.4020, -0.1433, -0.4672, -3.0769, 9.5483, 9.5043, -2.9127, -8.0081}; @ 1 to constrain any parameter to its starting value, else 0. @ CZVAR = 0; @ If CZVAR = 1, create a vector of zeros and ones of the same @ @ dimension as B identifying which parameters are to vary, else make @ @ BACTIVE = { 0 }. (column vector) @ BACTIVE = ones(8,1) | zeros(1,1) | ones(1,1); @ 1 to constrain any of the error components to equal each other, @ @ else 0. (integer) @ CEVAR = 0; @ If CEVAR=1, create a matrix of ones and zeros to constrain the @ @ necessary random parameters to be equal, else make EQMAT=0. (matrix)@ EQMAT = { 0 }; @ Identify the optimization routine: 1 Paul Ruud's Newton-Raph @ @ 2 for Gauss's maxlik @ OPTIM = 1; @ Specify the optimization algorithm/hessian calculation.(integer) @ METHOD = 2; @ 1 if OPTIM=1 and want to use modified BHHH, else 0. @ MIM = 1; @ If OPTIM=1, set STEP to the step length the maximization routine @ @ should initially try. @ STEP = 1; @ 1 to check the above inputs for conformability and consistency, @ @ else 0. @ VERBOSE = 1; @======================================================================@ @@@ You should not need to change anything below this line @@@ @======================================================================@ @ Create a series of needed globals @ NERR = NNC + NLC; @ Number of rnd draws @ NCOEF = NFC + NNC + NLC; @ Number of means @ NEVAR = rows(B); @ Length of beta @ @ Check inputs if VERBOSE=1 @ if ((VERBOSE /= 1) and (VERBOSE /= 0)); print "VERBOSE must 0 or 1."; print "Program terminated."; stop; endif; if VERBOSE; pcheck; else; print "The inputs have not been checked since VERBOSE=0."; print; endif; @ Rescale the variables. @ if RESCALE; j = rows(RSCLMAT); i = 1; if VERBOSE; if (RSCLMAT[.,1] > NVAR); print "RSCLMAT identifies a variable that is not in the data set."; print "Program terminated."; stop; endif; print "Rescaling Data:"; print " Variable Mult. Factor"; do while i <= j; RSCLMAT[i,1] RSCLMAT[i,2]; XMAT[.,(RSCLMAT[i,1]-1)*NALT+1:RSCLMAT[i,1]*NALT] = XMAT[.,(RSCLMAT[i,1]-1)*NALT+1:RSCLMAT[i,1]*NALT] * RSCLMAT[i,2]; i = i + 1; endo; print " "; else; do while i <= j; XMAT[.,(RSCLMAT[i,1]-1)*NALT+1:RSCLMAT[i,1]*NALT] = XMAT[.,(RSCLMAT[i,1]-1)*NALT+1:RSCLMAT[i,1]*NALT] * RSCLMAT[i,2]; i = i + 1; endo; endif; endif; @ Print out starting values if VERBOSE = 1. @ if VERBOSE; if NFC; print "The model has" NFC " fixed coefficients,"; print "for variables" IDC[1:NFC,1]; else; print "The model has no fixed coefficients."; endif; if NNC; print "The model has" NNC " normal coefficients,"; print "for variables" IDC[NFC+1:NFC+NNC,1]; else; print "The model has no normal coefficients."; endif; if NLC; print "The model has" NLC " fixed coefficients,"; print "for variables" IDC[NFC+NNC+1:NCOEF,1]; else; print "The model has no log-normal coefficients."; endif; print; if CZVAR; print "Starting values:" ; print B; print ; print "Parameters that are estimated (1) or held at starting value (0):"; print BACTIVE; print ; endif; if (CZVAR == 0); print "Starting values:"; print B; print; print "All parameters are estimated; none is held at its starting value."; print; endif; endif; @ Create new B and BACTIVE with reduced lengths if @ @ user specifies equality constraints. @ if CEVAR; if CZVAR; BACTIVE = EQMAT * BACTIVE; BACTIVE = BACTIVE .> 0; endif; B = (EQMAT * B) ./ (EQMAT * ones(NEVAR,1)); if VERBOSE and CZVAR; print "Some parameters are constrained to be the same."; print "Starting values of constrained parameters:"; print B; print ; print "Constrained parameters that are estimated (1) or held at starting value (0):"; print BACTIVE; print; endif; if VERBOSE and (CZVAR == 0); print "Some parameters are constrained to be the same."; print "Starting values of constrained parameters:"; print B; print; print "All constrained parameters are estimated; none is held at its starting value."; print; endif; endif; @ Describing random terms. @ if VERBOSE; print "Random error terms are based on:"; print "Seed: " SEED1; print "Repetitions: " NREP; print; print "Number of random draws for each observation for each repetitions: " (NERR); print; endif; @ Create dependent variable permutation matrix for the gradient @ YPERM=zeros(NOBS,NALT); i = 1; do while i<=NOBS; YPERM[i,YVEC[i,1]] = 1; i = i + 1; endo; @ Initialize the STEP to twice its value for Paul Ruud's routine. @ if (OPTIM==1); STEP = STEP*2; endif; @ Set up and do maximization @ if OPTIM == 1; {beta,ihesh} = domax(&ll,&gr,B,BACTIVE); Mix_b = beta; save Mix_b; Mix_cov = ihesh; save Mix_cov; print; print "Remember: Signs of standard deviations are irrelevant."; print "Interpret them as being positive."; print; endif; if OPTIM == 2; library maxlik,pgraph; #include maxlik.ext; maxset; _max_GradTol = EPS; _max_GradProc = &gr; _max_MaxIters = NITER; _max_Algorithm = METHOD; if DOWT; /* __weight = WEIGHT;*/ endif; if CZVAR; _max_Active = BACTIVE; endif; {beta,f,g,cov,ret} = maxlik(XMAT,0,&ll,B); call maxprt(beta,f,g,cov,ret); Mix_b = beta; save Mix_b; Mix_cov = cov; save Mix_cov; print; print "Remember: Signs of standard deviations are irrelevant."; print "Interpret them as being positive."; print; if (CZVAR == 0); print "gradient(hessian-inverse)gradient is:" ((g/_max_FinalHess)'g); else; g = selif(g,BACTIVE); print "gradient(hessian-inverse)gradient is:" ((g/_max_FinalHess)'g); endif; if PRTHESS; print "diagonal of hessian:" ((diag(_max_FinalHess))'); print; endif; if DOCH; ihesh = inv(_max_FinalHess); endif; endif; if NLC; mo = exp(beta[NFC+NNC+1:NCOEF,1]); me = exp(beta[NFC+NNC+1:NCOEF,1]+(beta[NCOEF+NNC+1:NEVAR,1]^2)/2); s = me.*( (exp(beta[NCOEF+NNC+1:NEVAR,1]^2) - 1).^(.5) ); print "Implied Means and Variances for Log-Normals"; print "--------------------------------------------------------------------------"; format /m1 /rds 12, 6; print " Coefficients Mode Mean Standard Dev"; print (seqa(1,1,NLC)~mo~me~s); print; endif; if DOCH; format /m1 /rds 12,8; L = zeros(NNC,NNC); if CZVAR; bt = beta.*BACTIVE; else; bt = beta; endif; j = NFC+NNC; k=1; do while k < NNC; L[k,.] = bt[j+1:j+k,1]' ~ zeros(1,NNC-k); j = j + k; k = k + 1; endo; L[k,.] = bt[j+1:j+k,1]'; V = L*L'; print "Variance-Covariance Matrix"; print "--------------------------------------------------------------------------"; print V; print; {Y,D} = stevar(L,ihesh[NFC+NNC+1:rows(ihesh),NFC+NNC+1:rows(ihesh)]); format /m1 /rds 6,3; print "Standard Errors"; print Y; print; print "t-statistics"; if CZVAR; print V./D; else; print V./Y; endif; print; format /m1 /rdt 16,8; endif; @ THIS IS THE END OF THE MAIN PROGRAM. @ @ PROCS par2coef, ll, gr, stevar, gpnr, expand, domax, pcheck follow.@ /* TURN THE NON-FIXED ELEMENTS OF THE PARAMETER VECTOR INTO AN NREPxq MATRIX OF COEFFICIENTS, WHERE q IS THE NUMBER OF NON-FIXED COEFFICIENTS TO BE ESTIMATED. */ proc (2) = par2coef(b); local coef, err, k, j; err = rndn(NREP, NERR); @ Seed must be set by calling proc @ if DOCH; coef = zeros(NREP,NERR); j = NNC; k = 1; do while k<=NERR; coef[.,k] = b[k,1]+sumc( (b[j+1:j+k,1]'.*err[.,1:k])' ); j = j + k; k = k + 1; endo; else; coef = b[1:NERR,1]'+b[NERR+1:NERR+NERR,1]'.*err; if NLC; coef[.,NNC+1:NERR] = exp(coef[.,NNC+1:NERR]); endif; endif; retp(coef, err); endp; /* LOG-LIKELIHOOD FUNCTION */ proc ll(b,x); @ Relies on the globals: NALT, NOBS, NP, NREP, SEED1, CENSOR @ @ NFC, NNC, NLC, IDC, NERR, NCOEF, NEVAR @ @ CZVAR, BACTIVE, CEVAR, EQMAT, XMAT @ local coef, c, k, n, t, km, kmm, rd, seed2; local v, err, ev, p0, p00; if CEVAR; b = EQMAT' * b; @ Expand b to its original size. @ endif; v = zeros(NOBS,NALT); @ Argument to logit formula @ p0 = zeros(NP,1); @ Simulated probability @ if CENSOR; c = (IDCENSOR-1)*NALT; @ Location of censor variable @ endif; k = 1; do while k <= NFC; @ Adds variables with fixed coefficients @ km = (IDC[k,1]-1)*NALT; v = v + b[k] .* XMAT[.,(km+1):(km+NALT)]; k = k+1; endo; seed2 = SEED1; rndseed seed2; rd = 0; n = 1; do while n <= NP; if NERR; {coef,err} = par2coef(b[NFC+1:NEVAR,1]); endif; p00 = ones(NREP,1); t=1; do while (t<=TIMES[n]); kmm = rd + t; ev = v[kmm,.]; k = 1; do while k <= NERR; @ Normal and Log-normal @ km = (IDC[NFC+k,1]-1)*NALT; ev = ev + coef[.,k] .* XMAT[kmm,(km+1):(km+NALT)]; k = k+1; endo; ev = exp(ev); if CENSOR; p00 = p00.*ev[.,YVEC[kmm,1]]./sumc((ev .* XMAT[kmm,(c+1):(c+NALT)])'); else; p00 = p00.*ev[.,YVEC[kmm,1]]./(sumc(ev')); endif; t = t+1; endo; rd = rd + TIMES[n]; p0[n,1] = meanc(p00); n = n + 1; endo; retp(ln(p0)); endp; /* GRADIENT PROCEDURE */ proc gr(b,x); @ Relies on the globals: NALT, NOBS, NP, NREP, SEED1, CENSOR @ @ NFC, NNC, NLC, IDC, NERR, NCOEF, NEVAR @ @ CZVAR, BACTIVE, CEVAR, EQMAT, XMAT, YPERM @ local coef, c, err, g, j, k, n, t, km, kmm, rd, seed2, y; local denom, der, v, ev, p0, p00, p1; if CEVAR; b = EQMAT' * b; @ Expand b to its original size. @ endif; if CENSOR; c = (IDCENSOR-1)*NALT; @ Location of censor variable @ endif; v = zeros(NOBS,NALT); @ Argument to logit formula @ p0 = zeros(NP,1); @ Simulated probability @ der = zeros(NP,NEVAR); @ Jacobian matrix @ k = 1; do while k <= NFC; @ Fixed coefficients @ km = (IDC[k,1]-1)*NALT; v = v + b[k] .* XMAT[.,(km+1):(km+NALT)]; k = k+1; endo; seed2 = SEED1; rndseed seed2; rd = 0; n = 1; do while n <= NP; if NERR; {coef,err} = par2coef(b[NFC+1:NEVAR,1]); endif; p00 = ones(NREP,1); g = zeros(NREP,NEVAR); t=1; do while (t<=TIMES[n]); kmm = rd + t; ev = v[kmm,.]; k = 1; do while k <= NERR; @ Normal and Log-normal @ km = (IDC[k+NFC,1]-1)*NALT; ev = ev + coef[.,k] .* XMAT[kmm,(km+1):(km+NALT)]; k = k+1; endo; ev = exp(ev); if CENSOR; denom = sumc((ev .* XMAT[kmm,(c+1):(c+NALT)])'); p00 = p00 .* ev[.,YVEC[kmm,1]]./denom; p1 = (ev.*XMAT[kmm,(c+1):(c+NALT)])./denom; else; denom = sumc(ev'); p00 = p00 .* ev[.,YVEC[kmm,1]]./denom; p1 = ev./denom; endif; y = YPERM[kmm,.]-p1; k = 1; do while k<=NCOEF; km = (IDC[k,1]-1)*NALT; g[.,k] = g[.,k] + sumc((XMAT[kmm,km+1:km+NALT].*y)'); k = k + 1; endo; t = t+1; endo; if DOCH; j = NFC+NNC; k = 1; do while k <= NNC; g[.,j+1:j+k] = g[.,NFC+k] .* err[.,1:k]; j = j + k; k = k + 1; endo; elseif NERR; if NLC; @ correction term for log-normal means @ g[.,NFC+NNC+1:NCOEF] = g[.,NFC+NNC+1:NCOEF] .* coef[.,NNC+1:NERR]; endif; @ errors for normals and log normals @ g[.,NCOEF+1:NEVAR] = g[.,NFC+1:NCOEF] .* err; endif; p0[n,1] = meanc(p00); der[n,.] = (meanc(p00.*g))'; rd = rd + TIMES[n]; n = n + 1; endo; if CEVAR; der = (der./p0) * EQMAT'; else; der = der./p0; endif; retp(der); endp; /* TURN THE STANDARD ERRORS OF THE CHOLESKY FACTOR INTO THE STANDARD ERRORS FOR THE ASSOCIATED VARIANCE-COVARIANCE MATRIX. */ proc (2) = stevar(L,ihesh); local CZD, CZihesh, cloc, rloc, D, k, r, s, seq; k = NNC; @ rows of L; number of variables in V @ s = k*k; @ rows(D); rows(vec(L)); rows(vec(V)) @ r = sumc(seqa(1,1,k)); @ cols of D; cols/rows of ihesh; @ D = zeros(s,r); @ derivative of L wrt L[j,i] @ cloc = 1; @ column location D @ j = 1; do while j<=k; i = 1; do while i<=j; rloc = (i-1)*k; seq = i; do while seq<=k; if (seq == j); D[rloc+1:rloc+k,cloc] = L[.,i]; D[rloc+j,cloc] = 2*D[rloc+j,cloc]; rloc = rloc+k; else; rloc = rloc+j-1; D[rloc+1,cloc] = L[seq,i]; rloc = rloc+1+(k-j); endif; seq = seq + 1; endo; cloc=cloc+1; i=i+1; endo; j = j+1; endo; if CZVAR; CZD = D .* BACTIVE[NFC+NNC+1:NEVAR]'; CZihesh = BACTIVE[NFC+NNC+1:NEVAR]' .* ihesh .* BACTIVE[NFC+NNC+1:NEVAR]; CZD = sqrt(diag(CZD*CZihesh*CZD')); D = sqrt(diag(D*ihesh*D')); @ standard errors @ retp(reshape(CZD,k,k), reshape(D,k,k)); else; D = sqrt(diag(D*ihesh*D')); @ standard errors @ retp(reshape(D,k,k), D); endif; endp; /* GRADIENT FOR PAUL RUUD'S ROUTINE WHEN USING NEWTON-RAPHSON*/ /* USE WHEN OPTIM == 1 AND METHOD == 2 */ proc gpnr(b); @ Relies on global: XMAT @ local grad; if DOWT; grad = WEIGHT .* gr(b,XMAT); else; grad = gr(b,XMAT); endif; retp(sumc(grad)); endp; /* EXPANDS THE DIRECTION VECTOR; ALLOWS PARAMETERS TO STAY AT STARTING */ /* VALUES; HELPER PROCEDURE FOR &domax */ proc expand( x, e ); local i,j; i = 0; j = 0; do while i < rows(e); i = i + 1; if e[i]; j = j + 1; e[i] = x[j]; endif; endo; if j/=rows(x); "Error in expand."; stop; endif; retp( e ); endp; /* MAXIMIZATION ROUTINE COURTESY OF PAUL RUUD */ proc (2) = domax( &f, &g, b, bactive ); @ Relies on the globals: CZVAR, EPS, METHOD, NITER, NOBS, @ @ PRTHESS, XMAT, NP @ local f:proc, g:proc; local direc, first, grad, sgrad, hesh, ihesh, lambda; local nb, printer, repeat, rhesh, step1, wtsq; local _biter, _f0, _f1, _fold, _tol; _tol = 1; _biter = 0; nb = seqa(1,1,rows(b)); _f0 = f( b, XMAT ); if DOWT; _f0 = sumc(WEIGHT.*_f0); wtsq = WEIGHT .^ (.5); else; _f0 = sumc(_f0); endif; format /m1 /rdt 16,8; print; print; print; do while (_tol > EPS or _tol < 0) and (_biter < NITER); _biter = _biter + 1; print "=========================================================================="; print " Iteration: " _biter; print " Function: " _f0; grad = g( b, XMAT ); if (METHOD == 1); if DOWT; grad = wtsq.*grad; hesh = grad'grad; grad = wtsq.*grad; else; hesh = grad'grad; endif; else; if DOWT; grad = WEIGHT .* grad; endif; hesh = -gradp( &gpnr, b ); @ WEIGHT done in &gpnr @ endif; sgrad = sumc(grad); @ Select only the variables that we want to maximize over @ if CZVAR; rhesh = hesh; sgrad = selif( sgrad, bactive ); hesh = selif( hesh, bactive ); hesh = selif( hesh', bactive ); endif; if (det(hesh)==0); print "Singular Hessian!"; print "Program terminated."; stop; else; ihesh = inv(hesh); direc = ihesh * sgrad; endif; _tol = direc'sgrad; if CZVAR; direc = expand( direc, bactive); endif; print " Tolerance: " _tol; print "--------------------------------------------------------------------------"; if PRTHESS; if CEVAR and CZVAR; printer = expand(sgrad./NP,bactive)~expand(diag(hesh),bactive); printer = nb~b~EQMAT'*printer[.,1]~EQMAT'printer[.,2]; elseif CEVAR; printer = nb~b~(EQMAT'*sgrad./NP)~(EQMAT'*diag(hesh)); elseif CZVAR; printer = nb~b~expand(sgrad./NP,bactive)~expand(diag(hesh),bactive); else; printer = nb~b~(sgrad./NP)~(diag(hesh)); endif; print " Coefficients Rel. Grad. Hessian"; else; if CEVAR and CZVAR; printer = expand(sgrad./NP,bactive); printer = nb~b~EQMAT'*printer[.,1]; elseif CEVAR; printer = nb~b~(EQMAT'*sgrad./NP); elseif CZVAR; printer = nb~b~expand(sgrad./NP,bactive); else; printer = nb~b~(sgrad./NP); endif; print " Coefficients Rel. Grad."; endif; print printer; if (_tol >= 0) and (_tol < EPS); break; elseif _tol < 0; direc = -direc; endif; step1 = STEP; lambda = .5; repeat = 1; first = 1; _f1 = _f0; steploop: step1 = step1 * lambda; _fold = _f1; if DOWT; _f1 = sumc(WEIGHT .* f(b+step1*direc,XMAT)); else; _f1 = sumc( f( b+step1*direc, XMAT ) ); endif; print "--------------------------------------------------------------------------"; print " Step: " step1; print " Function: " _f1; if repeat; print " Change: " _f1-_f0; else; print " Change: " _f1-_fold; endif; if MIM; if (step1 < 1.e-5); print "Failed to find increase."; retp(b); elseif (_f1 <= _f0) and (repeat); first = 0; goto steploop; elseif (_f1 > _fold) and (first); lambda = 2; repeat = 0; goto steploop; endif; else; if (step1 < 1.e-5); print "Failed to find increase."; retp(b); elseif (_f1 <= _f0); goto steploop; endif; endif; if (repeat); b = b + step1*direc; _f0 = _f1; else; b = b + .5 * step1 * direc; _f0 = _fold; endif; endo; print "=========================================================================="; print; format /m1 /rdt 1,8; if (_tol< EPS); print "Converged with tolerance: " _tol; print "Function value: " _f0; else; print "Stopped with tolerance: " _tol; print "Function value: " _f1; endif; print; lambda = eig(hesh); if lambda>0; print "Function is concave at stopping point."; else; print "WARNING: Function is not concave at stopping point."; endif; print; if DOWT; if CZVAR; grad = grad'grad; grad = selif(grad, BACTIVE); grad = selif(grad', BACTIVE); ihesh = ihesh*(grad)*ihesh; else; ihesh = ihesh*(grad'grad)*ihesh; endif; endif; if CEVAR and CZVAR; printer = expand(sqrt(diag(ihesh)), bactive); printer = nb~b~EQMAT'*printer; elseif CEVAR; printer = nb~b~(EQMAT'*sqrt(diag(ihesh))); elseif CZVAR; printer = nb~b~expand(sqrt(diag(ihesh)),bactive); else; printer = nb~b~sqrt((diag(ihesh))); endif; format /m1 /rdt 16,8; print " Parameters Estimates Standard Errors"; print "--------------------------------------------------------------------------"; print printer; if (CZVAR == 0); rhesh = hesh; endif; retp( b, inv(rhesh) ); endp; /* This proc checks the inputs if VERBOSE=1 */ proc (0) = pcheck; local i,j; @ Checking XMAT @ if (rows(XMAT) /= NOBS); print "XMAT has" rows(XMAT) "rows"; print "but it should have NOBS=" NOBS "rows."; print "Program terminated."; stop; elseif(cols(XMAT) /= (NVAR*NALT)); print "XMAT has" cols(XMAT) "columns"; print "but it should have NVAR*NALT= " (NVAR*NALT) "columns."; print "Program terminated."; stop; else; print "XMAT has:"; print "Rows: " NOBS; print "Cols: " NVAR*NALT ; print "Containing" NVAR "variables for " NALT " alternatives."; print; endif; @ Check dependent variable @ if (rows(YVEC) /= NOBS); print "YVEC has" rows(YVEC) "rows"; print "but it should have NOBS=" NOBS "rows."; print "Program terminated."; stop; endif; if (cols(YVEC) /= 1); print "YVEC should have only one column."; print "Program terminated."; stop; endif; @ Check TIMES @ if (rows(TIMES) /= NP); print "TIMES has" rows(TIMES) "rows"; print "but it should have NP=" NP "rows."; print "Program terminated."; stop; endif; if (cols(TIMES) /= 1); print "TIMES should have only one column."; print "Program terminated."; stop; endif; if (sumc(TIMES) /= NOBS); print "sumc of TIMES =" sumc(TIMES) "but it should equal NOBS."; print "Program terminated."; stop; endif; @ Check WEIGHT @ if (DOWT /=0) and (DOWT /=1); print "DOWT must be 0 or 1."; print "Program terminated."; stop; endif; if DOWT; if (rows(WEIGHT) /= NP); print "WEIGHT has" rows(WEIGHT) "rows"; print "but it should have NP=" NP "rows."; print "Program terminated."; stop; endif; if (cols(WEIGHT) /= 1); print "WEIGHT should have only one column."; print "Program terminated."; stop; endif; endif; @ Check NFC, NNC, NLC @ if ( (NFC+NNC+NLC) > NVAR ); print "NFC + NNC + NLC is greater than NVAR"; print "Program terminated."; stop; endif; @ Checking coefficients @ if (rows(IDC) /= NFC+NNC+NLC); print "IDC has" rows(IDC) "rows when it should have NFC+NNC+NLC=" NFC+NNC+NLC "rows."; print "Program terminated."; stop; endif; if (cols(IDC) /= 1); print "Commas are needed between the elements of IDC."; print "Program terminated."; stop; endif; if (sumc(IDC.>NVAR)>0); print "IDC identifies a variable that is not in the data set."; print "All elements of IDFC should be <= NVAR, which is " NVAR; print "Program terminated."; stop; endif; @ Check CENSOR @ if ((CENSOR /= 1) and (CENSOR /= 0)); print "CENSOR must be 0 or 1."; print "Program terminated."; stop; endif; if CENSOR; if IDCENSOR > NVAR; print "The censoring variable cannot be located"; print "since you set CENSOR larger than NVAR."; print "Program terminated."; stop; endif; i = (IDCENSOR-1)*NALT; j = 1; do while j<=NOBS; if (XMAT[j,i+YVEC[j,1]]==0); print "Your censoring variable eliminates the chosen alternative"; print "for observation " j; print "Program terminated."; stop; endif; j = j + 1; endo; endif; @ Check PRTHESS and RESCALE @ if ((PRTHESS /= 1) and (PRTHESS /= 0)); print "PRTHESS must be 0 or 1."; print "Program terminated."; stop; endif; if ((RESCALE /= 1) and (RESCALE /= 0)); print "RESCALE must be 0 or 1."; print "Program terminated."; stop; endif; @ Check DOCH @ if ((DOCH /= 1) and (DOCH /= 0)); print "DOCH must be 0 or 1."; print "Program terminated."; stop; endif; if DOCH; if (NLC /= 0); print "DOCH==1 and NLC is not zero."; print "The full var-cov matrix can be estimated with"; print "only normals and fixed coefficients."; print "Program terminated."; stop; endif; endif; @ Check B @ if (DOCH == 1); if (rows(B) /= NFC+NNC+sumc(seqa(1,1,NNC))); print "Starting values B has " rows(B) "rows"; print "when it should have NFC+NNC+(1+2+...+NNC)= " NFC+NNC+sumc(seqa(1,1,NNC)) "rows."; print "Program terminated."; stop; endif; print "Estimating the full var-cov matrix"; print; endif; if (DOCH==0) and (rows(B) /= NFC+2*NNC+2*NLC); print "Starting values B has " rows(B) "rows"; print "when it should have NFC+2*NNC+2*NLC= " NFC+2*NNC+2*NLC "rows."; print "Program terminated."; stop; endif; if (cols(B) /= 1); print "Commas needed between the elements of B."; print "Program terminated."; stop; endif; @ Check CZVAR @ if ((CZVAR /= 1) and (CZVAR /= 0)); print "CZVAR must be 0 or 1."; print "Program terminated."; stop; endif; if CZVAR; if (rows(BACTIVE) /= NEVAR); print "BACTIVE has " rows(BACTIVE) "rows"; print "when it should have NFC+2*NNC+2*NLC= " NEVAR "rows."; print "Program terminated."; stop; endif; if (cols(BACTIVE) /= 1); print "Commas needed between the elements of BACTIVE."; print "Program terminated."; stop; endif; endif; @ Check CEVAR @ if ((CEVAR /= 1) and (CEVAR /= 0)); print "CEVAR must be 0 or 1."; print "Program terminated."; stop; endif; if CEVAR; if (cols(EQMAT) /= NEVAR); print "EQMAT has " cols(EQMAT) " columns"; print "when it should have NFC+2*NNC+2*NLC=" NEVAR " columns."; print "Program terminated."; stop; endif; if (rows(EQMAT)>=NEVAR); print "EQMAT has " rows(EQMAT) " rows"; print "when it should have strictly less than NEVAR=" NEVAR; print "rows."; print "Program terminated."; stop; endif; endif; @ Checking NREP @ if (NREP <= 0); print "Error in NREP: must be positive."; print "Program terminated."; stop; endif; @ Check SEED1 @ if (SEED1<=0); print "SEED1 =" SEED1 "must be a positive integer."; print "Program terminated."; stop; endif; @ Check SV @ @ if (SV /= 0) and (SV /=1) and (SV /=2); print "SV=" SV "and it must be either 0, 1, or 2"; print "Program terminated."; stop; endif; if SV; if (cols(SEED1) /= 1); print "SEEDS must be a column vector."; print "Put commas between each element."; print "Program terminated."; stop; endif; endif; @ @ Check METHOD and OPTIM @ if (METHOD < 1); print "METHOD must be 1-6."; print "Program terminated."; stop; endif; if (OPTIM /= 1) and (OPTIM /= 2); print "OPTIM must be 1 or 2."; print "Program terminated."; stop; endif; if ((OPTIM == 1) and (METHOD > 2)); print "Method " METHOD " is not an option with OPTIM = 1."; print "Program terminted."; stop; endif; if ((OPTIM == 2) and (METHOD > 6)); print "Method " METHOD " is not an option with OPTIM = 2."; print "Program terminated."; stop; endif; @ Check MIM @ if ((MIM /= 1) and (MIM /= 0)); print "MIM must be 0 or 1."; print "Program terminated."; stop; endif; print "The inputs pass all the checks and seem fine."; print; retp; endp; closeall; end;