﻿

# 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