libname sasdata "C:\Users\pcname\Desktop\mysas";
Proc iml;
P=3 ; n=50 ; nruns=50; B= {1, 1, 1}; sum1=0;
X1= {0.0396637
0.7041423
0.5965423
3.7637131
1.9483883
0.6426752
4.5774973
2.8115605
1.7396808
-2.328225
-1.555311
1.9304524
-0.438061
-2.100597
4.1315962
-2.043585
1.7933469
-2.752522
2.3937606
1.6584021
-3.464741
-0.948967
-0.700051
0.9235539
1.0683674
0.1133843
-0.584628
-3.467575
2.480151
-0.23564
-0.449791
3.8318656
-1.594308
3.2146945
-1.972683
-2.409909
-1.222361
-1.850948
-2.457552
-2.786037
0.3688551
-0.745163
0.661344
0.9002574
-4.056465
-0.616641
-0.433525
0.893546
0.5791517
0.1619015};
x2={2.2098305
2.4159765
1.4874399
4.4362109
2.2857215
1.3777612
3.3380025
2.7249503
2.0605257
-4.323939
0.8696841
-0.296454
-1.572753
-4.024938
3.9134578
0.4531205
4.9395914
-6.149045
1.7525422
0.7839797
-5.230953
-6.2818
0.4636905
-2.611297
0.7944336
-3.29826
1.6261934
-3.250607
-0.745343
1.3766689
-3.611049
1.9873417
-1.146432
3.7802815
-0.053812
-4.514641
0.7478499
2.7953156
-1.689206
-4.303781
3.6997174
-0.822928
-1.584738
5.2773745
-4.61168
4.6989522
-4.669677
1.065509
-0.162652
3.5672716};
x3 = {2.5610138
2.4568404
3.3860373
4.2793987
2.1869485
3.9195027
5.1350544
1.2754415
4.0789311
-5.643419
4.2400767
-0.864723
-1.759528
-0.609483
3.9374213
6.3893745
2.8841345
-8.918494
-4.130928
-1.670308
-3.661193
-6.988873
1.7714041
2.1092845
4.00296
-5.230733
-2.968761
0.0599031
2.1487263
-1.713545
-3.31429
3.96951
-2.560647
3.0266561
-1.841822
-5.536233
-0.169626
2.4713955
-4.17196
-4.926937
1.6750868
1.7033197
-3.928748
4.4535524
-2.799376
1.2370943
-3.769872
2.9396264
-1.806647
5.5282356};
Xt1 = t(x1);
xt2=t(x2);
xt3=t(x3);
xt=xt1 || xt2||xt3;
Do j=1 to nruns;
E= normal(repeat(-1,n));
Y=xt* B+E;
my = y[:,];
m1 = xt1[:,];
m2 = xt2[:,];
m3 = xt3[:,];
ny =(y-my)/sqrt(t(y-my)*(y-my));
nX1 =(xt1-m1)/sqrt (t(xt1-m1)*(xt1-m1));
nX2 =(xt2-m2)/sqrt(t(xt2-m2)*(xt2-m2));
nX3 =(xt3-m3)/sqrt(t(xt3-m3)*(xt3-m3));
nx= nx1||nx2||nx3;
covar=xt`*xt/(nrow(xt)-1);
gdiag=diag(covar);
gg=sqrt(inv(gdiag));
R=gg*covar*gg;
val = eigval(R);
vec = eigvec(R);
call eigen (lambda, v, R);
vt = t(v);
beta=inv(nx`*nx)*nx`*ny;
beta2=beta##2;
g = vt * beta;
lenght = t(beta)*beta;
start CorrCoef(x, y) ;
xStd = standard(x) ;
yStd = standard(y) ;
df = nrow(x) - 1;
return ( xStd` * yStd / df ) ;
finish;
r1 = CorrCoef(ny, nx1) ;
r2 = CorrCoef(ny, nx2) ;
r3 = CorrCoef(ny, nx3) ;
rxy= r1// r2// r3;
thrd = t(rxy) * v ;
sixth = t(rxy) * v;
uhat = ny - nx*beta;
RSS=ssq(uhat);
estvar = RSS/(nrow(ny)-p);
col1=vt[,1];
col2=vt[,2];
col3=vt[,3];
reset storage = sasdata.mlibrary;
store estvar lambda col1 col2 col3 g;
nybar = sum(ny)/ nrow(ny);
TSS = ny`*ny-nrow(ny)*nybar**2;
ESS= TSS-RSS;
MSX= ESS/ (p-1);
Elkh1 = (estvar/beta2[1])+1/lambda[1];
Elkh2 = (estvar/beta2[2])+1/lambda[2];
Elkh3 = (estvar/beta2[3])+1/lambda[3];
Elkham= Elkh1// Elkh2// Elkh3;
Newx1=nx2||nx3;
Be1=inv(t(newx1)*newx1)*t(newx1)*nx1;
Hat1 = nx1-newx1*be1;
RS1=ssq(hat1);
Nx1bar = sum(nx1)/ nrow(nx1);
TS1 = nx1`*nx1-nrow(nx1)#nx1bar##2;
R21=1-RS1/TS1;
Vif1=1/1-r21;
Newx2=nx1||nx3;
Be2=inv(newx2`*newx2)*newx2`*nx2;
Hat2 = nx2-newx2*be2;
RS2=ssq(hat2);
Nx2bar = sum(nx2)/ nrow(nx2);
TS2 = nx2`*nx2-nrow(nx2)#nx2bar##2;
R22=1-RS2/TS2;
Vif2=1/1-r22;
Newx3=nx1||nx2;
Be3=inv(newx3`*newx3)*newx3`*nx3;
Hat3 = nx3-newx3*be3;
RS3=ssq(hat3);
Nx3bar = sum(nx3)/ nrow(nx3);
TS3 = nx3`*nx3-nrow(nx3)#nx3bar##2;
R23=1-RS3/TS3;
Vif3=1/1-r23;
Vif=vif1//vif2//vif3;
start fun(K);
sumf=0;
do i = 1 to 3;
reset storage = sasdata.mlibrary;
load estvar lambda g;
sumf = sumf + (estvar*lambda)+ (g* k)**2 /(lambda+ k)**2;
end;
return (sumf);
finish fun;
start con(k);
c=j(4,1,0);
sumc1=0;
do i = 1 to 3;
reset storage = sasdata.mlibrary;
load lambda col1 ;
sumc1 = sumc1 + (col1**2 *lambda/(lambda+k)**2);
end;
c[1]= 10-sumc1;
sumc2=0;
do i= 1 to 3;
reset storage = sasdata.mlibrary;
load lambda col2 ;
sumc2 = sumc2 + (col2**2 *lambda/(lambda+k)**2);
end;
c[2]= 10-sumc2;
sumc3=0;
do i= 1 to 3;
reset storage = sasdata.mlibrary;
load lambda col3 ;
sumc3 = sumc3 + (col3**2 *lambda/(lambda+k)**2);
end;
c[3]= 10-sumc3;
c[4]=(max(lambda)/min(lambda))-(max(lambda)+k)/(min(lambda)+k);
return (c);
finish con;
k=j(1,1,0);
optn=j(1,11,.); optn[4]=1; optn[6]=2; optn[11]=0; optn[10]=4;
call nlpqn (rc, kres, "fun", k, optn) nlc="con";
mse1=0;
do i = 1 to 3;
mse1 = mse1 + (estvar*lambda)+ (g* k)**2 /(lambda+ k)**2;
end;
Sum1 = sum1+ MSE1;
End;
Sumbar1 = sum1/nruns;
Print sumbar1;
quit;
