new; cls;
label={"A","B","C","D","E","F","G"};
mu={17,12,4,3,6,9,21};
sig={8.5,6,2,1.5,3,4.5,10.5};
times=1000;

fn model(x)=x[1,.]+x[2,.]+x[3,.]+x[4,.]+x[5,.]+x[6,.]+x[7,.];        /*  model setting */
call costmodel2(label,mu,sig,times);



/*
**  contrib2.txt - Calculation of Contributions from each variable by Rank Correlation.
** (C) Copyright 2004 Yosuke Amijima. All Rights Reserved.
**
**  Reference:  David Vose
**              "Risk Analysis: A Quantitative Guide" 2nd Edition, Ch.15, John Wiley & Sons 2000.
**
**
**  Purpose:    Computes contributions from each variable by rank correlation coefficients
**              assuming N(mu,sig) for each variable.
**
**  Format:     call costmodel2(label,mu,sig,times);
**
**  Input:      label    n x 1 vector, labels;
**
**              mu       n x 1 vector, mu's of 1st,2nd,3rd,..., nth variable
**
**              sig      n x 1 vector, sig's of 1st,2nd,3rd,..., nth variable
**
**              times    scalar, number of simulations
**
**  Output:     normrho  n x 1 vector, normalized rank correlation coefficients(contributions)
**
**
**  Notice:     You need to set the model you are trying to evaluate before calling this procedure
**              by vector notations, i.e. you must use '.*' and './' for multiplication and division.
**
**              Normalized spread is calculated with signs which would be the precise way when there
**              are some negative rank correlation coefficients. If you don't like that, replace
**              sig1=sig.*((-1)^(rho.<0)) to sig1=sig; .
*/

proc costmodel2(label,mu,sig,times);
   local n,dist,i,sum,rho,sig1,spread,normrho,mask,fmt;
   n=rows(mu);
   dist=zeros(n,times);
   i=1;
   do while i<=n;
      dist[i,.]=normal2(mu[i],sig[i],times)';
      i=i+1;
   endo;
   sum=model(dist);
   sum=sum';
screen off;
   rho=zeros(n,1);
   i=1;
   do while i<=n;
      rho[i]=speaman(sum,dist[i,.]');
      i=i+1;
   endo;
screen on;
   sig1=sig.*((-1)^(rho.<0));
   spread=sig1/sumc(sig1);
   normrho=rho/sumc(rho);
/* display the results */
   print "                                                          Normalized";
   print "   Component        Mean         s.d.  Rank corr.     Spread   Rank corr.";
   mask={0 1 1 1 1 1};
   fmt="*.*lg"~12~6;
   call printfm(label~mu~sig~rho~spread*100~normrho*100,mask,fmt);
   format/rz 11,6;
   print "                     SUM:" sumc(sig) sumc(rho) sumc(spread*100) sumc(normrho*100);
   print;
   format /mb1 /ros 16,8;
/* graph */
   library pgraph;
   graphset;
   pqgwin auto;
   ylabel("% Normalized Rank Correlation Coefficients");
   _pdate="";
   _protate=1;
   bar(0,normrho*100);   
   retp(normrho);
endp;

proc speaman(y,x);
   local n,u,v,rho;
   n=rows(x);
   u=ranksame(x);
   v=ranksame(y);
   rho=1-6*sumc((u-v)^2)/(n*(n^2-1));      /* (3.25) */
   print/rz x~y~u~v~(u-v)^2;
   print;
   print/rz "  n=" n;
   print/rz "rho=" rho;
   print;
   retp(rho);
endp;

proc speaman2(y,x);
   local n,u,v,SSuv,SSuu,SSvv,rho;
   n=rows(x);
   u=ranksame(x);
   v=ranksame(y);
   SSuv=sumc((u-meanc(u)).*(v-meanc(v)));
   SSuu=sumc((u-meanc(u))^2);
   SSvv=sumc((v-meanc(v))^2);
   rho=SSuv/sqrt(SSuu*SSvv);               /* (3.26) */
   print/rz x~y~u~v~(u-v)^2;
   print;
   print/rz "  n=" n;
   print/rz "rho=" rho;
   print;
   retp(rho);
endp;

proc ranksame(x);
   local n,y,i,j,k,index;
   n=rows(x);
   y=zeros(n,1);
   i=1;
   do while i<=n;
      index=sumc(x.==minc(x));
      k=0;
      j=0;
      do while j<=index-1;
         k=k+(i+j);
         j=j+1;
      endo;
      k=k/index;
      y=y+k*(x.==minc(x));
      x=x+(x.==minc(x))*1e256;
      i=i+index;
   endo;
   retp(y);
endp;




proc normal(mu,sig,n);
   local x;
   x=rndu(n,1);
/* transform from U(0,1) to Normal(mu,sig) */
   x=cdfni(x);
   x=mu+sig*x;
   retp(x);
endp;

proc normal2(mu,sig,n);
   local d,U,i,x,index;
/* divide between 0 and 1 into n segments and sample from each segment */
   d=1/n;
   U=zeros(n,1);
   i=1;
   do while i<=n;
      U[i]=(d*i-d*(i-1))*rndu(1,1)+d*(i-1);
      i=i+1;
   endo;
/* randomize the order */
   x=zeros(n,1);
   i=1;
   do while i<=n;
      index=ceil((n-(i-1))*rndu(1,1));
      x[i]=U[index];
      U=delif(U,U.==U[index]);
      i=i+1;
   endo;
/* transform from U(0,1) to Normal(mu,sig) */
   x=cdfni(x);
   x=mu+sig*x;
   retp(x);
endp;