new; cls; /* The simulation of multivariate empirical distribution based on stratified sampling */ let data[200,3]= -0.66 -0.34 0.52 -0.24 0.43 0.43 0.59 1.07 0.67 0.22 0.32 0.89 2.32 2.63 1.99 0.32 0.68 0.16 -0.98 -0.96 -0.7 -1.94 -1.68 -1.64 -0.67 -1.15 -1.03 0.9 0.58 0.64 -1.08 -0.95 -1.3 0.12 0.61 -0.26 -1.45 -1.18 -1.34 1.48 1.17 2.03 -0.41 -0.63 -0.14 -0.13 0.18 -0.05 1.43 0.74 1.19 2.24 2.44 1.92 0.04 -0.19 0.61 0.45 0.5 0.33 -0.31 -0.09 -0.91 1.05 0.78 0.39 0.82 -0.32 0.64 -0.08 0.4 0.23 1.15 0.55 1.08 -0.02 -0.29 0.57 -0.66 -0.6 -1.01 0.73 0.77 0.92 -0.14 0.08 -0.16 0.07 0.51 -0.05 1.49 1.66 1.43 -0.65 -0.39 -0.66 0.22 0.7 0.78 -0.34 -0.87 -0.58 0.65 0.41 0.42 0.16 0.39 0.17 -0.37 0.88 0.5 0.73 1.02 0.97 1.29 1.79 1.11 0.65 0.48 0.65 0.7 0.8 0.81 -1.43 -1.07 -1.29 -0.12 0.27 -0.09 -0.61 -0.28 0.16 -0.99 -0.94 -1.18 -2.09 -1.96 -1.79 0.29 0.23 0.82 1.53 1.67 1.82 0.48 -1.07 -0.84 1.15 0.91 0.96 0.78 0.54 0.94 -0.19 -0.77 -0.85 1.51 2.01 1.99 0.38 0.36 0.38 -0.13 -0.8 0.01 0.3 -0.32 0.13 0.5 0.15 0.35 1.68 1.24 1.13 0.44 0.19 0.31 0.67 1.42 0.37 1.12 0.6 0.7 0.01 0.3 0.03 0.66 0.68 0.41 0.76 0.56 1.09 1.19 1.59 1.22 -1.26 -0.87 -1.11 -0.36 -0.73 0.36 0.25 0.53 -0.4 0.28 0.28 0.62 -0.49 -0.57 -0.48 -0.44 -0.38 0.02 -1.06 -0.66 -1 0.18 0.04 0.41 -0.49 0.23 0.41 0.42 0.2 0.32 -0.44 -0.02 -0.28 0.24 0.22 0.41 -1.08 -1.38 -1.48 0.35 0.24 0.63 0.25 0.54 0.84 -1.38 -1.52 -0.76 -1.88 -1.57 -1.43 0.77 1.01 1.64 -0.7 -0.43 -1.07 -0.15 0.25 0.14 -0.86 -0.45 -0.66 -1.3 -0.91 -0.99 0.39 0.29 0.2 0.72 0.14 0.77 -0.53 -0.72 -0.42 -1.42 -1.8 -0.73 0.22 0.59 0.51 0.83 1.65 0.31 -1.89 -1.89 -2.52 -1.06 -1.36 -1.44 -0.51 0.14 -0.05 -1.78 -1.54 -1.34 -1.87 -1.39 -1.75 -1.31 -1.54 -1.6 0.64 0.39 0.48 -0.76 -0.66 -0.75 -0.57 -0.67 -0.26 -0.69 -0.26 -0.76 -0.35 -0.19 -0.19 -0.28 0.1 0.19 -1.22 -1.12 -0.87 0.2 -0.14 0.03 1.82 1.38 1.27 -1.71 -1.11 -1.41 0.2 0.63 0.43 -0.51 -1.09 -1.21 0.33 0.28 0.59 -0.4 -0.43 0.1 -0.19 -0.4 -0.84 0.6 0.92 1.14 -0.17 -0.24 -0.11 1.13 1.54 1.4 -1.21 -0.83 -0.59 -1.26 -1.05 -1.2 -1.76 -1.26 -1.47 0.2 0.21 0.62 -0.37 -1.15 -0.62 -0.58 -0.92 -0.47 0.92 0.16 -0.12 -0.97 -1.17 -0.7 -0.6 -0.8 -1.11 -0.87 -0.54 -0.81 -0.34 -0.22 -0.72 0.98 1.54 1.46 1.58 1.41 1.2 -0.62 -1.31 -0.67 0.02 0.25 0.59 0.03 0.08 0.05 -0.8 -0.3 -0.67 -0.27 0.05 -0.14 -0.79 -0.67 -0.4 0.3 0.27 0.25 0.97 1.36 1.3 0.18 -0.22 0.59 0.07 -0.29 -0.52 0.91 1.03 0.8 -0.73 -0.8 -0.54 0.4 -0.12 0.28 -0.41 -0.65 -0.05 -1.3 -0.3 -0.28 0.36 0.31 0.36 0.33 0.25 0.52 0.72 -0.05 0.26 1.36 1.38 1.76 0.82 0.13 0.25 -0.88 -0.63 -0.62 0.27 0.2 0.27 -1.42 -0.86 -0.24 -0.5 -1.1 -0.25 -1 -0.82 -0.7 1.65 1.21 2.39 -0.03 -0.7 -0.64 0.85 0.16 0.21 0.12 0.29 0.54 1.76 1.8 2.3 0.88 0.22 0.6 -2.74 -2.46 -2.14 0.86 0.36 1.29 -0.37 -1.16 -0.92 1.98 1.3 1.08 0.28 0.48 0.77 0 0.14 0.14 -0.33 -1.12 0.26 0.06 -0.3 0.38 -0.67 -1.42 -1.17 -0.08 0.16 0.24 -0.47 -0.93 -0.4 -0.91 -0.56 -0.75 0.86 1.43 0.63 -0.95 -0.9 -1.13 0.83 1 1.1 -0.79 -1.64 -0.83 -2.12 -1.89 -1.17 -0.73 -0.55 -0.59 -1.27 -0.59 -0.21 -2.69 -2.22 -2.05 1.33 1.86 1.43 0.14 0.08 0.54 -0.5 0.07 -0.53 -0.8 -0.73 -0.8 1.44 1.43 1.18 1.16 -0.11 0.65 0.41 1.75 1.46 0.15 -0.16 -0.05 -0.97 -1.61 -1.63 -0.09 0.27 0.14 -0.19 -0.37 -0.07 1.52 1.42 1.91 0.93 1.31 0.99 1.27 1.49 1.1 0.11 -0.46 0.76 -0.88 -0.93 -1.33 -0.4 0.4 0.51 1.62 1.6 0.86 0.42 -0.22 0.18 ; Y=empdistextsim2s(data); library pgraph; graphset; pqgwin auto; _plctrl=-1; title("original distribution"); xlabel("X"); ylabel("Y"); zlabel("Z"); xyz(data[.,1],data[.,2],data[.,3]); title("simulated distribution"); xyz(Y[.,1],Y[.,2],Y[.,3]); /* ** mempextsims.txt - Multivariate Empirical Monte Carlo Sampler considering ties ** based on stratified sampling with tails extended. ** (C) Copyright 2007 Yosuke Amijima. All Rights Reserved. ** ** Purpose: Samples multivariate empirical Monte Carlo through ** the multivariate empirical Copula generated from data. ** ** Format: Y=empdistextsim2s(data); ** ** Input: data matrix, n x k matrix of original data ** ** Output: Y matrix, n x k matrix of simulated ** ** ** Notice: The range of each series of the simulated data are going to become ** the same as that of the original series on average. The tails are ** greater than or less than those of the original series. This ** simulation holds the rank patterns of data so the Kendall's tau ** and Spearman's rho are preserved. We should not apply this program ** to the data with some boundaries. */ proc empdistextsim2s(data); local n,k,U,Y,i,Ui,yi,index; n=rows(data); k=cols(data); U=empsim2s(data); Y=zeros(n,k); i=1; do while i<=k; Ui=sortc(U[.,i],1); Ui=delif(Ui,(Ui[1:n-1].==Ui[2:n])|0); yi=rndempexti2(data[.,i],Ui); yi=sortc(yi,1); index=ranksame0(data[.,i]); Y[.,i]=yi[index]; i=i+1; endo; retp(Y); endp; proc empsim2s(data); local n,k,index,sindex,nn,i,j,U,Ui,Ui1; n=rows(data); k=cols(data); U=zeros(n,k); i=1; do while i<=k; index=ranksame0(data[.,i]); sindex=sortc(index,1); nn=maxc(index); Ui=U01ss(n); Ui=sortc(Ui,1); Ui1=zeros(nn,1); j=1; do while j<=nn; Ui1[j]=sumc(Ui.*(sindex.==j))/sumc(sindex.==j); j=j+1; endo; U[.,i]=Ui1[index]; i=i+1; endo; retp(U); endp; proc U01ss(nn); /* not randomizing the order */ local d,i,U; /* divide between 0 and 1 into n segments and sample from each segment */ d=1/nn; U=zeros(nn,1); i=1; do while i<=nn; U[i]=(d*i-d*(i-1))*rndu(1,1)+d*(i-1); i=i+1; endo; /* check if the first element is zero or not */ do while U[1]==0; U[1]=d*rndu(1,1); endo; retp(U); endp; proc ranksame0(x); local n,index,x0,nn,y,i; n=rows(x); x0=sortc(x,1); index=(x0[1:n-1].==x0[2:n]); x0=delif(x0,index|0); nn=rows(x0); y=zeros(n,1); i=1; do while i<=n; y[i]=(x[i].==x0)'seqa(1,1,nn); i=i+1; endo; retp(y); endp; proc rndempexti2(x,u); local n,i,y,yy,p,x0,n0,xm,j,index; n=rows(x); x=sortc(x,1); i=ranksameS(x); i=sortc(i,1); i=delif(i,(i[1:n-1].==i[2:n])|0); p=0|(i/n); x0=delif(x,(x[1:n-1].==x[2:n])|0); n0=rows(x0); xm=(x0[1]-(x0[2]-x0[1])/2)|((x0[1:n0-1]+x0[2:n0])/2)|(x0[n0]+(x0[n0]-x0[n0-1])/2); yy=zeros(rows(u),1); j=1; do while j<=rows(u); if u[j]==1; y=xm[n0+1]; else; index=sumc(u[j].>=p); /* calculate the index number=1,2,3,...,n0 */ y=xm[index]+(xm[index+1]-xm[index])*rndu(1,1); endif; yy[j]=y; j=j+1; endo; retp(yy); endp; proc ranksameS(x); local n,x0,index,y,i,s; n=rows(x); x0=sortc(x,1); index=rankindx(x,1); y=zeros(n,1); i=1; do while i<=n; s=maxc((x.==x0[i]).*index); y=y+s.*(x.==x0[i]); i=s+1; endo; retp(y); endp;