option nonotes; options validvarname = v6; ***This option is needed if you are running the program in SAS8. Since I wrote the program using SAS6.12, I used the CALIS variable intercep (instead of intercept). Without this options statement, SAS8 will not recognize intercep. Note, if you are running the program using SAS 6.12, this option can be deleted; data param; reliab = .75; gam0= 2 ; gam1= 1 ; gam2= 1 ; gam3= 2; b10 = 1 ; b20 = 2 ; b30 = 3 ; b40 = 4 ; b50 = 5 ; b60 = 6 ; b11 = .7 ; b21 = .3 ; b32 = .4 ; b42 = .8 ; b53 = .5 ; b63 = .4 ; ps1= .5; ps2=.5 ; ps3=.5; ps4=.5 ; ps5=.5 ; ps6=.5 ; ps7=.5 ; ps8=.5 ; ps9=.5 ; mf1 = .5; mf2 = -.5; ph1= 1 ; ph12=.5; ph2 = 1; zeta= .3 ; ps1 = sqrt(b11*b11*ph1*(1-reliab)/reliab); ps2 = sqrt(b21*b21*ph1*(1-reliab)/reliab); ps3 = sqrt(b32*b32*ph2*(1-reliab)/reliab); ps4 = sqrt(b42*b42*ph2*(1-reliab)/reliab); ps7 = sqrt(ph1*(1-reliab)/reliab); ps8 = sqrt(ph2*(1-reliab)/reliab); varf1f2 = (ph1+mf1*mf1)*(ph2+mf2*mf2) + (ph12+mf1*mf2)**2 - 2*mf1*mf1*mf2*mf2; cov1f1f2= mf1*ph12+ph1*mf2; cov2f1f2= mf2*ph12+ph2*mf1; phf3 = gam1*gam1*ph1 + gam2*gam2*ph2 + gam3*gam3*varf1f2 + 2*gam1*gam3*cov1f1f2 + 2*gam2*gam3*cov2f1f2 + 2*ph12 + zeta*zeta; phf3 =7.9; ps5 = sqrt(b53*b53*phf3*(1-reliab)/reliab); ps6 = sqrt(b63*b63*phf3*(1-reliab)/reliab); ps9 = sqrt(phf3*(1-reliab)/reliab); run; proc iml; use param; read all var{ph1 ph2 ph12} into temp; phi=j(2,2,0); out=j(1,3,0); phi[1,1]=temp[1,1]; phi[2,2]=temp[1,2]; phi[1,2]=temp[1,3]; phi[2,1]=temp[1,3]; halfphi= half(phi); out[1,1]= halfphi[1,1];out[1,2]= halfphi[1,2]; out[1,3]= halfphi[2,2]; create param1 from out[colname = {hph11 hph12 hph22}]; append from out; close param1; quit; data param2; merge param param1; run; proc datasets nolist; delete saveparm save; run; %macro docalis(sam,obs); data start; do i=1 to %eval(&sam*&obs); *xf1=rannor(125); xf1=ranuni(125); xf1=sqrt(12)*(xf1-.5); *xf2=rannor(125); xf2=ranuni(125); xf2=sqrt(12)*(xf2-.5); %do j=3 %to 12; x&j=rannor(125); %end; output; end; data a; set start; if _n_=1 then set param2; f1= mf1 + hph11*xf1; f2= mf2 + hph12*xf1 + hph22*xf2; f3= gam0 + gam1*f1 + gam2*f2 + gam3*f1*f2 + zeta*x3; Y1 = b10+ b11*f1 + ps1*x4; Y2 = b20+ b21*f1 + ps2*x5; Y3 = b30+ b32*f2 + ps3*x6; Y4 = b40+ b42*f2 + ps4*x7; Y5 = b50+ b53*f3 + ps5*x8; Y6 = b60+ b63*f3 + ps6*x9; Y7 = f1 + ps7*x10; Y8 = f2 + ps8*x11; Y9 = f3 + ps9*x12; output; keep Y1 Y2 Y3 Y4 Y5 Y6 Y7 Y8 Y9 f1 f2; run; data a; set a; %do p=1 %to &sam; if %eval(1+(&p-1)*&obs)<=_n_<=%eval(&p*&obs) then iter=&p; %end; proc means data=a noprint; by iter; var y1-y9; output out=outmean mean=mean1-mean9; run; data a; set a; merge a outmean; by iter; run; PROC CALIS METHOD=ML omethod=quanew uCOV augment DATA=a MAXITER=5000 maxfunc=2000 outram=modelone outest=covtest noprint; by iter; var y1-y9; LINEQS Y1 = m1 intercep + B11 F1 + E1, Y2 = m2 intercep +B21 F1 + E2, Y3 = m3 intercep +B32 F2 + E3, Y4 = m4 intercep +B42 F2 + E4, Y5 = m5 intercep +B53 F3 + E5, Y6 = m6 intercep +B63 F3 + E6, Y7 = m7 intercep + F1 + E7, Y8 = m8 intercep + F2 + E8, Y9 = m9 intercep + F3 + E9; STD E1-E9 = TH1-TH9, F1-F3=PH1-PH3; COV F1-F3=PH4-PH6; BOUNDS 0<=PH1-PH3, 0<=TH1-TH9; parameters b10, b20, b30, b40, b50, b60; m1 = m7*b11 + b10; m2 = m7*b21 + b20; m3 = m8*b32 + b30; m4 = m8*b42 + b40; m5 = m9*b53 + b50; m6 = m9*b63 + b60; data saveparm(keep = iter b10 b20 b30 b40 b50 b60 B11 B21 B32 B42 B53 B63 TH1-TH9 PH1-PH6 seb10 seb20 seb30 seb40 seb50 seb60 seb11 seb21 seb32 seb42 seb53 seb63 seth1 seth2 seth3 seth4 seth5 seth6 seth7 seth8 seth9); set modelone; if _NAME_="B10" then do;B10 = _estim_;seB10 = _stderr_;end; if _NAME_="B20" then do;B20 = _estim_;seB20 = _stderr_;end; if _NAME_="B30" then do;B30 = _estim_;seB30 = _stderr_;end; if _NAME_="B40" then do;B40 = _estim_;seB40 = _stderr_;end; if _NAME_="B50" then do;B50 = _estim_;seB50 = _stderr_;end; if _NAME_="B60" then do;B60 = _estim_;seB60 = _stderr_;end; if _NAME_="B11" then do;B11 = _estim_;seB11 = _stderr_;end; if _NAME_="B21" then do;B21 = _estim_;seB21 = _stderr_;end; if _NAME_="B32" then do;B32 = _estim_;seB32 = _stderr_;end; if _NAME_="B42" then do;B42 = _estim_;seB42 = _stderr_;end; if _NAME_="B53" then do;B53 = _estim_;seB53 = _stderr_;end; if _NAME_="B63" then do;B63 = _estim_;seB63 = _stderr_;end; if _NAME_="TH1" then do;TH1 = _estim_;seTH1 = _stderr_;end; if _NAME_="TH2" then do;TH2 = _estim_;seTH2 = _stderr_;end; if _NAME_="TH3" then do;TH3 = _estim_;seTH3 = _stderr_;end; if _NAME_="TH4" then do;TH4 = _estim_;seTH4 = _stderr_;end; if _NAME_="TH5" then do;TH5 = _estim_;seTH5 = _stderr_;end; if _NAME_="TH6" then do;TH6 = _estim_;seTH6 = _stderr_;end; if _NAME_="TH7" then do;TH7 = _estim_;seTH7 = _stderr_;end; if _NAME_="TH8" then do;TH8 = _estim_;seTH8 = _stderr_;end; if _NAME_="TH9" then do;TH9 = _estim_;seTH9 = _stderr_;end; if _NAME_="PH1" then do;PH1 = _estim_;end; if _NAME_="PH2" then do;PH2 = _estim_;end; if _NAME_="PH3" then do;PH3 = _estim_;end; if _NAME_="PH4" then do;PH4 = _estim_;end; if _NAME_="PH5" then do;PH5 = _estim_;end; if _NAME_="PH6" then do;PH6 = _estim_;end; /*find out what ph4,5,6 represent see below */ run; proc means noprint; by iter; var b10 b20 b30 b40 b50 b60 b11 b21 b32 b42 b53 b63 th1 th2 th3 th4 th5 th6 th7 th8 th9 ph1 ph2 ph3 ph4 ph5 ph6 seb10 seb20 seb30 seb40 seb50 seb60 seb11 seb21 seb32 seb42 seb53 seb63 seth1 seth2 seth3 seth4 seth5 seth6 seth7 seth8 seth9; output out= new(keep= iter b10 b20 b30 b40 b50 b60 b11 b21 b32 b42 b53 b63 th1 th2 th3 th4 th5 th6 th7 th8 th9 ph1 ph2 ph3 ph4 ph5 ph6 seb10 seb20 seb30 seb40 seb50 seb60 seb11 seb21 seb32 seb42 seb53 seb63 seth1 seth2 seth3 seth4 seth5 seth6 seth7 seth8 seth9) mean = b10 b20 b30 b40 b50 b60 b11 b21 b32 b42 b53 b63 th1 th2 th3 th4 th5 th6 th7 th8 th9 ph1 ph2 ph3 ph4 ph5 ph6 seb10 seb20 seb30 seb40 seb50 seb60 seb11 seb21 seb32 seb42 seb53 seb63 seth1 seth2 seth3 seth4 seth5 seth6 seth7 seth8 seth9; run; data savecov(keep = iter b10 b20 b30 b40 b50 b60 b11 b21 b32 b42 b53 b63 th1 th2 th3 th4 th5 th6 th7 th8 th9); set covtest; if (_TYPE_ = "COV") then do; if (_NAME_ = "B10") then output; if (_NAME_ = "B20") then output; if (_NAME_ = "B30") then output; if (_NAME_ = "B40") then output; if (_NAME_ = "B50") then output; if (_NAME_ = "B60") then output; if (_NAME_ = "B11") then output; if (_NAME_ = "B21") then output; if (_NAME_ = "B32") then output; if (_NAME_ = "B42") then output; if (_NAME_ = "B53") then output; if (_NAME_ = "B63") then output; if (_NAME_ = "TH1") then output; if (_NAME_ = "TH2") then output; if (_NAME_ = "TH3") then output; if (_NAME_ = "TH4") then output; if (_NAME_ = "TH5") then output; if (_NAME_ = "TH6") then output; if (_NAME_ = "TH7") then output; if (_NAME_ = "TH8") then output; if (_NAME_ = "TH9") then output; end; data save; set new; merge new outmean; by iter; run; %mend docalis; %docalis(1000,500); **********This (1000,500) is for running the simulation 1000 times and using data with sample size = 500; data cov; set savecov; run; proc iml; use save; read all var{TH1 TH2 TH3 TH4 TH5 TH6 TH7 TH8 TH9} into thetac; read all var{B11 B21 B32 B42 B53 B63} into betac; read all var{PH1 PH2 PH3 PH4 PH5 PH6} into phic; read all var{mean1 mean2 mean3 mean4 mean5 mean6 mean7 mean8 mean9} into means; use param1; read point 1 var{hph11 hph12 hph22} into hph; use cov; read all var{b10 b20 b30 b40 b50 b60 b11 b21 b32 b42 b53 b63 th1 th2 th3 th4 th5 th6 th7 th8 th9} into covdata; itn=1000; n=500; **********This itn = 1000 and n = 500 is for running the simulation 1000 times and using data with sample size = 500; reliab = .75; gam0= 2 ; gam1= 1 ; gam2= 1 ; gam3= 2; b10 = 1 ; b20 = 2 ; b30 = 3 ; b40 = 4 ; b50 = 5 ; b60 = 6 ; b11 = .7 ; b21 = .3 ; b32 = .4 ; b42 = .8 ; b53 = .5 ; b63 = .4 ; ps1= .5; ps2=.5 ; ps3=.5; ps4=.5 ; ps5=.5 ; ps6=.5 ; ps7=.5 ; ps8=.5 ; ps9=.5 ; mf1 = .5; mf2 = -.5; ph1= 1 ; ph12=.5; ph2 = 1; zeta= .3 ; ps1 = sqrt(b11*b11*ph1*(1-reliab)/reliab); ps2 = sqrt(b21*b21*ph1*(1-reliab)/reliab); ps3 = sqrt(b32*b32*ph2*(1-reliab)/reliab); ps4 = sqrt(b42*b42*ph2*(1-reliab)/reliab); ps7 = sqrt(ph1*(1-reliab)/reliab); ps8 = sqrt(ph2*(1-reliab)/reliab); varf1f2 = (ph1+mf1*mf1)*(ph2+mf2*mf2) + (ph12+mf1*mf2)**2 - 2*mf1*mf1*mf2*mf2; cov1f1f2= mf1*ph12+ph1*mf2; cov2f1f2= mf2*ph12+ph2*mf1; phf3 = gam1*gam1*ph1 + gam2*gam2*ph2 + gam3*gam3*varf1f2 + 2*gam1*gam3*cov1f1f2 + 2*gam2*gam3*cov2f1f2 + 2*ph12 + zeta*zeta; phf3=7.9; ps5 = sqrt(b53*b53*phf3*(1-reliab)/reliab); ps6 = sqrt(b63*b63*phf3*(1-reliab)/reliab); ps9 = sqrt(phf3*(1-reliab)/reliab); do it =1 to itn; y1=j(n,1); y2=y1; y3=y1; y4=y1; y5=y1; y6=y1; y7=y1; y8=y1; y9=y1; int=y1; do i=1 to n; *x1=rannor(125);*x2=rannor(125); x1=ranuni(125); x1=sqrt(12)*(x1-.5); x2=ranuni(125); x2=sqrt(12)*(x2-.5); x3=rannor(125);x4=rannor(125); x5=rannor(125);x6=rannor(125);x7=rannor(125);x8=rannor(125); x9=rannor(125);x10=rannor(125);x11=rannor(125);x12=rannor(125); f1= mf1 + hph[1,1]*x1; f2= mf2 + hph[1,2]*x1 + hph[1,3]*x2; f3= gam0 + gam1*f1 + gam2*f2 + gam3*f1*f2 + zeta*x3; Y1[i,] = b10+ b11*f1 + ps1*x4; Y2[i,] = b20+ b21*f1 + ps2*x5; Y3[i,] = b30+ b32*f2 + ps3*x6; Y4[i,] = b40+ b42*f2 + ps4*x7; Y5[i,] = b50+ b53*f3 + ps5*x8; Y6[i,] = b60+ b63*f3 + ps6*x9; Y7[i,] = f1 + ps7*x10; Y8[i,] = f2 + ps8*x11; Y9[i,] = f3 + ps9*x12; end; y=y1||y2||y3||y4||y5||y6||y7||y8||y9; psi= diag(thetac[it,]); beta=j(6,3,0); beta[1,1]=betac[it,1]; beta[2,1]=betac[it,2]; beta[3,2]=betac[it,3]; beta[4,2]=betac[it,4]; beta[5,3]=betac[it,5]; beta[6,3]=betac[it,6]; beta=beta//i(3); phi=j(3,3,0); phi[1,1]=phic[it,1]; phi[2,2]=phic[it,2]; phi[3,3]=phic[it,3]; phi[1,2]=phic[it,4]; phi[2,1]=phic[it,4]; phi[1,3]=phic[it,5]; phi[3,1]=phic[it,5]; phi[2,3]=phic[it,6]; phi[3,2]=phic[it,6]; ybar=means[it,1:6]; xbar=means[it,7:9]; covarian = j(21,21,0); do k=1 to 6; row= 21*(it-1)+15+k; covarian[k,]=covdata[row,]; end; do k=1 to 15; row1 = 6 + k; row2 = 21*(it-1)+k; covarian[row1,]=covdata[row2,]; end; beta0=j(9,1,0); hold= ybar`- beta[1:6,]*xbar`; beta0= hold//j(3,1,0); beta0=hold; bigbeta0=beta0@j(1,n); sigvv=psi[1:6,1:6]+(beta[1:6,]*psi[7:9,7:9]*(beta[1:6,])`); siguv=-psi[7:9,7:9]*(beta[1:6,])`; varrho=psi[7:9,7:9]-siguv*inv(sigvv)*(siguv)`; invsigvv=inv(sigvv); v=y[,1:6]`-bigbeta0-beta[1:6,]*y[,7:9]`; *fhat=inv((beta)`*inv(psi)*beta)*(beta)`*inv(psi)*(Y`-bigbeta0); fhat = (y[,7:9])` - siguv * invsigvv * v; k=ncol(beta); fhat12=fhat[1:2,]; f1f2=fhat[1,]#fhat[2,]; ff1=fhat[1,]#fhat[1,]; ff2=fhat[2,]#fhat[2,]; hh=j(n,1)||fhat12`||f1f2`; mhat=hh`*hh/n; fhat3= (fhat[3,])`; mhyhat=hh`*fhat3/n; sigmauu=j(4,4,0); sigmauu[2,2]= varrho[1,1]; sigmauu[2,4]=xbar[1,2]* varrho[1,1]; sigmauu[3,3]= varrho[2,2]; sigmauu[3,4]=xbar[1,1]* varrho[2,2]; sigmauu[4,2]=sigmauu[2,4]; sigmauu[4,3]=sigmauu[3,4]; *sigmauu[4,4]=(phi[1,1]+xbar[1,1]*xbar[1,1])*varrho[2,2]+(phi[2,2]+xbar[1,2]*xbar[1,2])*varrho[1,1]+varrho[1,1]*varrho[2,2]; sigmauu[4,4]=(sum(ff1)/n)*varrho[2,2]+(sum(ff2)/n)*varrho[1,1]-varrho[1,1]*varrho[2,2]; gamma=inv(mhat-sigmauu)*mhyhat; /*calculating the standard errors */ littlem = hh`# (fhat3`//fhat3`//fhat3`//fhat3`); subtract = j(4,n,0); do t=1 to n; sigmat =j(4,4,0); sigmat[2,2]= varrho[1,1]; sigmat[2,4]=fhat[2,t]* varrho[1,1]; sigmat[3,3]= varrho[2,2]; sigmat[3,4]=fhat[1,t]* varrho[2,2]; sigmat[4,2]=sigmat[2,4]; sigmat[4,3]=sigmat[3,4]; sigmat[4,4]=ff1[1,t]*varrho[2,2]+ff2[1,t]*varrho[1,1]-varrho[1,1]*varrho[2,2]; bigh=(hh[t,])`*hh[t,] - sigmat; subtract[,t]=bigh*gamma; end; dt=littlem - subtract; square1 = j(4,4,0); square2 = j(4,4,0); square1[2,2]=1; square1[2,4]=sum(fhat[2,])/n; square1[4,4]=sum(ff2)/n - varrho[2,2]; square1[4,2]=sum(fhat[2,])/n; square2[3,3]=1; square2[3,4]=sum(fhat[1,])/n; square2[4,4]=sum(ff1)/n - varrho[1,1]; square2[4,3]=sum(fhat[1,])/n; dldsig1 = square1*gamma; dldsig2 = square2*gamma; dsigth = i(3) - beta[1:6,]`*invsigvv*beta[1:6,]*psi[7:9,7:9] + siguv*invsigvv*beta[1:6,] + siguv*invsigvv*(beta[1:6,]*beta[1:6,]`)*invsigvv*siguv`; dsig1th7 = dsigth[1,1]; dsig2th7 = 0; dsig2th8 = dsigth[2,2]; dsig1th8 = 0; /*don't need dsig3. Also dsig1, dsig2 are zero w.r.t. th9 */ dsig1th1 = trace(siguv*invsigvv*diag(1||0||0||0||0||0)*invsigvv*siguv`); dsig2th1= 0; dsig1th2 = trace(siguv*invsigvv*diag(0||1||0||0||0||0)*invsigvv*siguv`); dsig2th2= 0; dsig2th3 = trace(siguv*invsigvv*diag(0||0||1||0||0||0)*invsigvv*siguv`); dsig1th3 =0; dsig2th4 = trace(siguv*invsigvv*diag(0||0||0||1||0||0)*invsigvv*siguv`); dsig1th4 =0; outside=j(3,6,0); outside[1,]=psi[7,7]*(1||0||0||0||0||0); dsg1b11 = trace(outside * invsigvv * siguv` + siguv*invsigvv*outside` + siguv*invsigvv*(beta[1:6,1]*(1||0||0||0||0||0) + (1//0//0//0//0//0)*beta[1:6,1]`)*psi[7,7]*invsigvv*siguv`); dsg2b11 = 0; outside=j(3,6,0); outside[1,]=psi[7,7]*(0||1||0||0||0||0); dsg1b21 = trace(outside * invsigvv * siguv` + siguv*invsigvv*outside` + siguv*invsigvv*(beta[1:6,1]*(0||1||0||0||0||0) + (0//1//0//0//0//0)*beta[1:6,1]`)*psi[7,7]*invsigvv*siguv`); dsg2b21 = 0; outside=j(3,6,0); outside[2,]=psi[8,8]*(0||0||1||0||0||0); dsg2b32 = trace(outside * invsigvv * siguv` + siguv*invsigvv*outside` + siguv*invsigvv*(beta[1:6,2]*(0||0||1||0||0||0) + (0//0//1//0//0//0)*beta[1:6,2]`)*psi[8,8]*invsigvv*siguv`); dsg1b32 = 0; outside=j(3,6,0); outside[2,]=psi[8,8]*(0||0||0||1||0||0); dsg2b42 = trace(outside * invsigvv * siguv` + siguv*invsigvv*outside` + siguv*invsigvv*(beta[1:6,2]*(0||0||0||1||0||0) + (0//0//0//1//0//0)*beta[1:6,2]`)*psi[8,8]*invsigvv*siguv`); dsg1b42 = 0; dsig1all=(j(1,6,0)||dsg1b11||dsg1b21||dsg1b32||dsg1b42||j(1,2,0)|| dsig1th1||dsig1th2||dsig1th3||dsig1th4||j(1,2,0)||dsig1th7||dsig1th8||0); dsig2all=(j(1,6,0)||dsg2b11||dsg2b21||dsg2b32||dsg2b42||j(1,2,0)|| dsig2th1||dsig2th2||dsig2th3||dsig2th4||j(1,2,0)||dsig2th7||dsig2th8||0); *dsig1all=j(1,21,0); *dsig2all=j(1,21,0); coeff = (i(12)||j(12,9,0))//dsig1all//dsig2all; covnew = coeff*covarian*coeff`; biggam = siguv*invsigvv; average = j(4,12,0); do t = 1 to n; part1a=j(4,4,0); part2a=j(4,4,0); part3a=j(4,1,0); part1b=j(4,4,0); part2b=j(4,4,0); part3b=j(4,1,0); part1c=j(4,1,0); part1a[1,2]=1; part1a[1,4]=fhat[2,t]; part1a[2,2]=2*fhat[1,t]; part1a[2,3]=fhat[2,t]; part1a[2,4]=2*f1f2[1,t]; part1a[3,4]=ff2[1,t]; part1a[4,4]=2*fhat[1,t]*ff2[1,t]; part1a[2,1]=1; part1a[4,1]=fhat[2,t]; part1a[3,2]=fhat[2,t]; part1a[4,2]=2*f1f2[1,t]; part1a[4,3]=ff2[1,t]; part2a[3,4]=varrho[2,2]; part2a[4,4]=2*fhat[1,t]*varrho[2,2]; part2a[4,3]=varrho[2,2]; inside1 = part1a -part2a; part3a[2,1]=1; part3a[4,1]=fhat[2,t]; dldfa = part3a*fhat[3,t] - inside1*gamma; part1b[1,3]=1; part1b[1,4]=fhat[1,t]; part1b[2,3]=fhat[1,t]; part1b[2,4]=ff1[1,t]; part1b[3,3]=2*fhat[2,t]; part1b[3,4]=2*f1f2[1,t]; part1b[4,4]=2*fhat[2,t]*ff1[1,t]; part1b[3,1]=1; part1b[4,1]=fhat[1,t]; part1b[3,2]=fhat[1,t]; part1b[4,2]=ff1[1,t]; part1b[4,3]=2*f1f2[1,t]; part2b[2,4]=varrho[1,1]; part2b[4,4]=2*fhat[2,t]*varrho[1,1]; part2b[4,2]=varrho[1,1]; inside2 = part1b -part2b; part3b[3,1]=1; part3b[4,1]=fhat[1,t]; dldfb = part3b*fhat[3,t] - inside2*gamma; part1c=hh[t,]`; dldf = dldfa||dldfb||part1c; b2t=biggam*diag(y[t,7]||y[t,7]||y[t,8]||y[t,8]||y[t,9]||y[t,9]); betat = biggam||b2t; dldfbeta=dldf*betat; average=average + dldfbeta#(1/n); end; climit=average||dldsig1||dldsig2; addition = climit*covnew*climit`; omegah = dt* dt` # (1/n) #(1/n) + addition; vhat = inv(mhat-sigmauu)*omegah*inv(mhat-sigmauu); segam0=sqrt(vhat[1,1]); segam1=sqrt(vhat[2,2]); segam2=sqrt(vhat[3,3]); segam3=sqrt(vhat[4,4]); myy=fhat3`*fhat3/n; cm=(myy||mhyhat`)//(mhyhat||mhat); rcm=(inv(root(cm)))`; sigma0=j(1,5,0)//(j(4,1,0)||sigmauu); msm=rcm*sigma0*rcm`; ll=eigval(msm); l=1/ll[1,]; lk=1><(l-(1/n)); a3=lk-3/n; bhatn3=inv(mhat-sigmauu#a3)*mhyhat; a8=lk-8/n; bhatn8=inv(mhat-sigmauu#a8)*mhyhat; a11=lk-1/n; bhatn1=inv(mhat-sigmauu#a11)*mhyhat; result = result//(gamma`||segam0||segam1||segam2||segam3||bhatn3`||bhatn8`||bhatn1`); end; create res from result; append from result; close res; quit; run; data relabel;set res; rename col1=gam0 col2=gam1 col3=gam2 col4=gam3 col5=segam0 col6=segam1 col7=segam2 col8=segam3 col9=bhatn30 col10=bhatn31 col11=bhatn32 col12=bhatn33 col13=bhatn80 col14=bhatn81 col15=bhatn82 col16=bhatn83 col17=bhatn10 col18=bhatn11 col19=bhatn12 col20=bhatn13; run; data finale; merge relabel save; drop mean1-mean6; run;