Original SAS code in Bharath and Shumway (2008 RFS) for Merton DD. See comment below for some minor issues.

Shumway.sas

/* This SAS program calculates the distance to default using the KMV-Merton model 
   with the iterated estimate of the volatility of firm value.  Many of the results 
   of Bharath and Shumway (2004) are generated by this program.  The program 
   requires the data described below, and it generates a permanent sas data file
   called ssd.kmv which contains distances to default.  The program calculates 
   monthly distances to default every year from 1980 to 1990 as it is currently written*/

/* This program requires two datasets: 

    ssd.comp,  which contains monthly observations taken from the quarterly compustat 
               data file for compustat items 45 and 51 (data45 and data51, respectively), 
               the crsp permno firm identifier, the calendar year and the calendar month 
               (calyr and calmnth, respectively) and the risk-free rate of return (r, 
               taken from some source other than compustat).

    ssd.dailycrsp, which has daily observations of crsp permno, current shares outstanding
               (shrout), the current date in sas format (date), and the price (prc). 
 */
 
libname ssd 'c:\kmv.dir';
data sampl;  set ssd.comp;
   cdt = 100*calyr + calmnth;
   f = 1000*(data45 + 0.5 * data51);
   if data45 < 0 or data51 < 0 then delete;
   drop data45 data51 calyr calmnth;

proc sort;  by permno cdt;

data ssd.kmv99; curdat = 0;

%macro itera(yyy,mmm);
  data one;
   set ssd.dailycrsp;
    if (100*(&yyy-1) + &mmm) <= (100*year(date) + month(date)) <= (100*&yyy + &mmm); 
    e = abs(prc)*shrout;
    cdt = 100*year(date) + month(date);
    keep permno date cdt e;

  proc sort;    by permno cdt;

  data one;
   merge one sampl;
    by permno cdt;
    if e ne . and f > 0;
    a = e + f;
    if a = . or permno = . or e = 0 then delete;

  proc sort;    by permno date;

* get volatility of total asset returns and look for large values;

   data one;   set one;   l1p = lag1(permno);   l1a = lag1(a);
   data one;   set one;   if l1p = permno then ra = log(a/l1a);
   proc means noprint data = one;   var ra f e;   by permno;   output out = bob;
   data bob1;   set bob;   if _stat_ = 'STD'; if _freq_ < 50 then delete; 
        va = sqrt(252)*RA;  if va < .01 then va = .01;    keep permno va;
   data bob2; set bob; if _stat_ = 'MEAN'; if f > 100000 and e > 100000 then largev = 1; 
        else largev = 0; keep largev permno;
     
  data one;   merge one bob1 bob2;   by permno;     if va = . then delete;
   if largev = 1 then do; f = f/10000; e = e/10000; a = a/10000; end;
  drop ra l1a l1p;

  data conv;   permno = 0;

*iteration;

 %do j = 1 %to 15;
   proc model noprint data = one;  endogenous a;  exogenous r f VA e;
      e = a*probnorm((log(a/f) + (r+va*va/2))/VA) - 
      f*exp(-r)*probnorm((log(a/f) + (r-va*va/2))/VA);  solve a/out=two;
   data two;   set two;   num = _n_;   keep a num;
   data one;   set one;   num = _n_;   drop a;  
   data two;   merge one two;   by num;  l1p = lag1(permno);   l1a = lag1(a);
   data two;   set two; if l1p = permno then ra = log(a/l1a);
   proc means noprint data = two;   var ra;   by permno;   output out = bob;
   data bar;   set bob; if _stat_ = 'MEAN';  mu = 252*ra;   keep permno mu;
   data bob;   set bob; if _stat_ = 'STD'; va1 = sqrt(252)*RA;
               if va1 < 0.01 then va1 = 0.01; keep permno va1;
   data one;   merge two bob bar; by permno; vdif = va1 - va; 
               if abs(vdif) < 0.001 and vdif ne . then conv = 1;
   data fin;   set one; if conv = 1;   assetvol = va1; proc sort; by permno descending date;  
   data fin;   set fin; if permno ne lag1(permno);   curdat = 100*&yyy + &mmm;  iter = &j;
   data conv;   merge conv fin;   by permno;  drop va ra l1p l1a conv cdt num;
   data one;   set one; if conv ne 1;   va = va1;   drop va1; 
 %end;

  data ssd.kmv;
   merge ssd.kmv conv;
    by curdat;
    if permno = 0 or curdat = 0 then delete;  drop va1;
    edf = probnorm(-((log(a/f) + (mu-(assetvol**2)/2))/assetvol));
    label edf = 'expected default frequency';
    label curdat = 'date in yyyymm format';
    label e = 'market equity';
    label iter = 'iterations required';
    label assetvol = 'volatility of a';
    label f = 'current debt + 0.5LTD';
    label vdif = 'assetvol - penultimate VA';
    label a = 'total firm value';
    label r = 'risk-free rate';
    label largev = 'one if assets, equity and f deflated';
    label mu = 'expected asset return';
%mend;

%macro bob;
    %do i = 1980 %to 2003;
        %do m = 1 %to 12;
             %itera(&i,&m);
        %end;
    %end;
%mend;
%bob;

There are two minor issues in this version of code:

1. The initial value of $\sigma_V$ is not set to $\sigma_E[E/(E+F)]$ as described in the paper.
2. It does not use the past year's data but the past month.
   • At line 36, cdt=100*year(date)+month(date) accidentally restricts the "past year" daily stock returns to the "past month" later. Note that at line 42-43 it merges by permno and cdt, where cdt refers to a certain year-month. We can pause the program after this data step to confirm that indeed there is only a month of data for each permno.
   • A correction is to change line 36 to cdt=100*&yyy.+&mmm.;.

For a self-contained complete code, see https://mingze-gao.com/posts/merton-dd