pmbrown

NearestCorr.sas

Nov 11th, 2020
30
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. /******************************************************************************
  2. Author, date: Rick Wicklin, NOV2012
  3. Macro name:   NearestCorr
  4. Description:  implement the Higham projection method to obtain the nearest
  5.               correlation matrix to a given symmetric correlation/covariance matrix
  6. Reference:    http://blogs.sas.com/content/iml/2012/11/28/
  7.               computing-the-nearest-correlation-matrix.html
  8. ******************************************************************************/
  9.  
  10. %macro NearestCorr();
  11. /* Project symmetric X onto S={positive semidefinite matrices}.
  12.    Replace any negative eigenvalues of X with zero */
  13. start ProjS(X);
  14.    call eigen(D, Q, X); /* note X = Q*D*Q` */
  15.    V = choose(D>.0001, D, .0001); /*pmbrown: this line changed as per comment
  16.                                    at bottom of sas blog*/
  17.    W = Q#sqrt(V`);      /* form Q*diag(V)*Q` */
  18.    return( W*W` );      /* W*W` = Q*diag(V)*Q` */
  19. finish;
  20.  
  21. /* project square X onto U={matrices with unit diagonal}.
  22.    Return X with the diagonal elements replaced by ones. */
  23. start ProjU(X);
  24.    n = nrow(X);
  25.    Y = X;
  26.    diagIdx = do(1, n*n, n+1);
  27.    Y[diagIdx] = 1;      /* set diagonal elements to 1 */
  28.    return ( Y );
  29. finish;
  30.  
  31. /* Helper function: the infinity norm is the max abs value of the row sums */
  32. start MatInfNorm(A);
  33.    return( max(abs(A[,+])) );
  34. finish;
  35.  
  36. /* Given a symmetric correlation matrix, A,
  37.    project A onto the space of positive semidefinite matrices.
  38.    The function uses the algorithm of Higham (2002) to return
  39.    the matrix X that is closest to A in the Frobenius norm. */
  40. start NearestCorr(A);
  41.    maxIter = 100; tol = 1e-8;  /* parameters...you can change these */
  42.    iter = 1;      maxd   = 1;  /* initial values */
  43.    Yold = A;  Xold = A;  dS = 0;
  44.  
  45.    do while( (iter <= maxIter) & (maxd > tol) );
  46.      R = Yold - dS; /* dS is Dykstra's correction */
  47.      X = ProjS(R);  /* project onto S={positive semidefinite} */
  48.      dS = X - R;
  49.      Y = ProjU(X);  /* project onto U={matrices with unit diagonal} */
  50.  
  51.      /* How much has X changed? (Eqn 4.1) */
  52.      dx = MatInfNorm(X-Xold) / MatInfNorm(X);
  53.      dy = MatInfNorm(Y-Yold) / MatInfNorm(Y);
  54.      dxy = MatInfNorm(Y - X) / MatInfNorm(Y);
  55.      maxd = max(dx,dy,dxy);
  56.  
  57.      iter = iter + 1;
  58.      Xold = X; Yold = Y; /* update matrices */
  59.    end;
  60.    return( X ); /* X is positive semidefinite */
  61. finish;
  62. %mend NearestCorr;
  63.  
  64. ***end********************************************************;
RAW Paste Data