new 200,50000; format /rd 8,4; outwidth 250; output file = simple.out reset; nmoments=66; /* Load in residual covs (cov2.sas => rescov1.dat */ load t1[] = rescov-male.dat; rescov = reshape(t1,nmoments,1); clear t1; /* Load in var-cov of residual covs */ load t2[] = vrescov-male.dat; vrescov = reshape(t2,nmoments,nmoments); clear t2; print "min distance model for residual covariances -- YOUNG MEN"; print "number of moments = " nmoments; print "vector of input moments and their sampling errors"; check1=seqa(1,1,nmoments)~rescov~sqrt(diag(vrescov)); print check1; clear check1; nparm=3; print "number of parameters = " nparm; b0={.04, .02, .8}; bnames={"parm1","parm2","parm3"}; check2=seqa(1,1,nparm)~b0; print "starting values"; print check2; clear check2; /* global parameters*/ _qn_RelGradTol= 5e-6; /*max relative gradient tolerance, def=1e-5 */ _qn_GradProc=0; /*use numerical gradient*/ _qn_MaxIters=500; /*max iterations*/ _qn_PrintIters=0; /*print iteration data 1=yes 0=no*/ _qn_ParNames=bnames; _qn_PrintResults=1; /*print results 1=y 0=no*/ print "Goodness of fit of Initial Inputs - Check of Coding"; test=fit(b0); print "test sum of squares= " test; f=bigf(b0); print "Jacobian = F matrix"; b1=f[.,1:3]; print b1; /* now fit model, print out stats, and calculate cov matrix */ {bex,vex,gex,retex}=QNewton(&ssq, b0); print "Goodness of fit of Final Estimates"; test=fit(bex); f=bigf(bex); h=f' * f; hinv=invpd(h); vbeta=hinv*f'*vrescov*f*hinv; se=sqrt(diag(vbeta)); check4=bex~se; print "beta and standard errors"; print check4; /* procedures = ltr, fit, ssq, bigf */ proc ltr(i,j); local even; even= ( i/2 - floor(i/2) )==0; retp ( even * ( (i-1)*(i-2)/2 + i-1 + j ) + (1-even) * ( i*(i-1)/2 + j ) ); endp; proc fit(b); local vomega, vu, rho, i,j,index,mfit,dev,actfit, se, t, ssr, rmse; vomega=b[1]; vu=b[2]; rho=b[3]; mfit=zeros(nmoments,1); i=1; index=1; do while i<=11; j=1; do while j