/****mweak nearly sing. design with Nese Yildiz*****/ /****for table 2, cp=10,m=15, rho=0.5, k=1-1/n, n=100, these can be changed to get other table values****/ ite=10000; "no of iterations:" ite; pb=zeros(ite,1);phhs=zeros(ite,1);wald=zeros(ite,1);cnum=zeros(ite,1);minei=zeros(ite,1);v2=zeros(ite,1); cls; iii=1; do while iii<=ite; n=100; beta0=0;/***tval***/ eta=rndn(n,1); v=rndn(n,1); r=0.5;/***corre between equation****/ m=15;/***no. of instruments****/ z=rndn(n,m); om=eye(m);k=1-1/n; om[1,2]=k; om[2,1]=k; z=z*chol(om);/***correlated instruments***/ e=(r*eta)+(sqrt(1-(r^2))*v); cp=10; pp=sqrt(cp/(((2*k)+m)*n))*ones(m,1);/****reduced form coeffs****/ /*** x is n by 1*****/ /*** z is n by m****/ x=z*pp+eta; y=x*beta0+e; /****************/ library optmum; d=zeros(5,1); h=zeros(5,1);hhs=zeros(5,1); bet0={-2,-1,0,1,2}; i=1; do while i<=5; __output=0; _opmiter=50; {b,f,g,retcode}=optmum(&cue,bet0[i]); h[i]=abs(f);/*** I also tried without abs , the results donot change, Abs is used since gauss sometimes gives negative nos, this is done to avoid that****/ d[i]=b; i=i+1; endo; j=minindc(h); pb[iii]=d[j];/***optimal b****/ zx=z.*x; ghat=meanc(zx); ze=z.*(y-x*pb[iii]); mze=meanc(ze); omega=ze'ze/n; /****OMEGA IN THE PAPER*****/ dlab=(z.*x)'(z.*x)/n; ei=eig(omega); cnum[iii]=maxc(ei)/minc(ei);/****CONDITION NUMBER*****/ minei[iii]=minc(ei);/****MIN EIGENVALUE****/ dzx=zx'ze/n; dhat=-meanc(zx)+(dzx*invpd(omega)*meanc(ze));/*** THIS IS ZHAT(BETAHAT) ON P.14 OF OUR ARTICLE****/ v1=n*dhat'invpd(omega)*dhat;/***** THIS IS THE MIDDLE TERM IN VAR COVAR: ZHAT' (INV OMEGA) ZHAT :SHOWN IN EQN 8 IN THE PAPER******/ hes0=ghat'invpd(omega)*ghat; hesa1=2*ghat'invpd(omega)*(-dzx-dzx')*invpd(omega)*mze; hesa2=mze'invpd(omega)*(-dzx-dzx')*invpd(omega)*(-dzx-dzx')*invpd(omega)*mze; hesa3=-mze'invpd(omega)*dlab'invpd(omega)*mze; hes1=n*(hes0+hesa1+hesa2+hesa3);/****ANALTYICAL HESSIAN TERM******/ v2[iii]=invpd(hes1)*v1*invpd(hes1);/*****VHAT IN ARTICLE IN WALD TEST , P.14-15******/ wald[iii]=(pb[iii]^2)/(v2[iii]); iii; iii=iii+1; endo; pb=sortc(pb,1); "median bias : " pb[.5*ite,1]; "Condition Number : " meanc(cnum); "Min eigvalue : " meanc(minei); v2=sortc(v2,1); "median std error : " ((v2[.5*ite,1]^0.5)); iqr=pb[.75*ite,1]-pb[.25*ite,1]; "interquartile range : " iqr; rej=(wald.>3.84); "size of the wald at 5% level : " meanc(rej); /***estimation***/ proc cue(bet); local g,om,of; g=z.*(y-x*bet); trap 1; om=invpd(g'g/n); trap 0; of=meanc(g)'om*meanc(g)/2; retp(of); endp;