A subroutine in SAS to simulate asset pricing paths

For matrix computation in SAS, SAS/IML is the choice. This module has its own syntax, functions and even plotting subsystem. Some statisticians used it to realize the algorithms beyond the reach of SAS’s procedures, for example, boosting [Ref. 1]. However, comparing with other popular matrix-based languages, such as R and Matlab, SAS/IML has no edge. SAS’s most valuable products are still its robust data step and statistical procedures. ‘Porting’ source codes from other languages into SAS has to rely on data step.

Asset prices can be estimated by Monte Carlo simulation. To generate a series of price-evolving paths with several most common parameters, Dr. Brandimarte codes a naive Matlab function to apply the standard Wiener process [Ref. 2]. With random seeds for normal distribution, multiple pricing mechanisms can be demonstrated and compared. In SAS 9.2, the workflow of simulation and visualization would be modularized as a data step subroutine. Later, such a pricing subroutine could be easily invoked under given circumstances.

It is said that SAS 9.3 is going to be released 2011Q3. Hope this time, the data step function compiler, Proc Fcmp, could be more dynamic and with more methods.

References:
1. Dmitrienko, Alex, Christy Chuang-Stein, and Ralph D’Agostino. Pharmaceutical Statistics Using SAS: A Practical Guide. SAS Publishing. 2007
2. Paolo Brandimarte. Numerical methods in finance. John Wiley & Sons. 2002.

/*******************READ ME*********************************************
* - A SUBROUTINE IN SAS TO SIMULATE ASSET PRICING PATHS -
*
* VERSION:     SAS 9.2(ts2m0), windows 64bit
* DATE:        17apr2011
* AUTHOR:      hchao8@gmail.com
*
****************END OF READ ME*****************************************/

****************(1) MODULE-BUILDING STEP******************;
******(1.1) COMPILE FUNCTION-ACCOMPANYING MACRO*************;
option mstored sasmstore = sasuser;
%macro AssetPath_macro() / store source;
   data _tmp1;
      set _tmp;
      day  = _n_ - 1;
   run;

   proc transpose data = _tmp1 out = _tmp2 ;
      by day;
      var path:;
   run;

   data _tmp2;
      set _tmp2;
      label _name_ = 'Simulated paths';
   run;

   ods html style = money;
   proc sgplot data = _tmp2;
      series x = day y = col1 / group = _name_;
      yaxis grid label = 'Asset price';
      xaxis grid label = 'Change by days';
   run;
   ods html close;
%mend;

******(1.2) COMPILE SUBROUTINE FOR ASSET PRICING*************;
proc fcmp outlib = sasuser.astpth.funcs;
   subroutine AssetPath(S0,  mu, sigma, T,  NSteps,  NRepl);
   /*************************************************************
   * SUBROUTINE:    AssetPath()
   * PARAMETERS:    S0    = the initial price
   *                mu    = the drift
   *                sigma = the volatility
   *                T     = the horizontal time
   *                NSteps= the number of time steps
   *                NRepl = the number of replications
   *************************************************************/
      array SPaths[1, 1] / nosymbols;
      array Path[1, 1] / nosymbols;
      call dynamic_array(SPaths, NRepl, NSteps + 1);
      call dynamic_array(Path, NSteps + 1, NRepl);
      call zeromatrix(SPaths);
      do _row = 1 to NRepl;
         SPaths[_row, 1] = S0;
      end;
      dt = T / NSteps;
      nudt = (mu - 0.5*sigma**2) * dt;
      sidt = sigma * sqrt(dt);
      do _row = 1 to NRepl;
         do _col = 1 to Nsteps;
            SPaths[_row, _col + 1] = SPaths[_row, _col] * 
               exp(nudt + sidt*rannor(0));
         end;
      end;   
      call transpose(SPaths, Path);
      rc1 = write_array('_tmp', Path);
      rc2 = run_macro('AssetPath_macro');
    endsub;
quit;
****************END OF STEP (1)******************;

****************(2) TESTING STEP******************;
option cmplib = (sasuser.astpth) mstored sasmstore = sasuser;
data _null_;
   S0 = 50; mu = 0.1; sigma = 0.3; T = 1; NSteps = 365; NRepl = 3;
   call AssetPath(S0,  mu, sigma, T,  NSteps,  NRepl);
run;
****************END OF STEP (2)******************;

****************END OF ALL CODING***************************************;

Predict unemployment rate for Election 2012 by SAS

Since recently President Obama announced that he is seeking reelection, the unemployment rate on November 2012 would decide the result. The Wall Street Journal averaged 54 economists’ predication and concluded that the number is going to be 7.7%. Apparently, those economists rely on the historical data to forecast the future, together with more or less their subjective judgment. However, the newly released March data is surprisingly good: 8.8%, which means that this predication number has to be adjusted downwardly to be below 7.7%. Then what is the real-time prediction of the unemployment rate for this ‘big’ time?

SAS has one of the finest time-series packages in the world: SAS/ETS which includes a few predictive procedures such as the ARIMA procedure and the FORECAST procedure[Ref. 2]. And the economic data is updated by Federal Reserve and well accessible on their website. To predict unemployment rate like a real professional is possible with a notebook computer and SAS. Of course SAS’s procedures have tons of methods and parameters to tune. To simply this problem, in the SAS macro below, I chose a conservative method and an aggressive one, to give a rough estimation about the unemployment range. Just like what the WSJ said, the trend matters. The predication will be more approaching to the real number as time goes forward. Right now, my prediction for the unemployment rate on November 2012 is from 7.1% to 7.4%.

References:
1."Jobless Rate at 2012 Presidential Vote Forecast at 7.7%, Highest Since Carter-Ford, but the Trend May Matter Most". The Wall Street Journal, 13Mar2011.
2.SAS/ETS 9.2 User Guide. SAS Publishing, 2008.

/*******************READ ME*********************************************
* -- PREDICT UNEMPLOYMENT FOR ELECTION 2012 LIKE A PRO --
*
* VERSION:     SAS 9.2(ts2m0), windows 64bit
* DATE:        09apr2011
* AUTHOR:      hchao8@gmail.com
*
****************END OF READ ME*****************************************/

****************(1) MODULE-BUILDING STEP******************;
%macro unemrate(predtime = );
   /***********************************************************
   *  MACRO:      unemrate()
   *  GOAL:       use time series based on latest FED data to
   *              predict unemployment rate in US and plot
   *  PARAMETERS: predtime = the time when unemployement rate 
   *                         is to be predicted
   *
   ***********************************************************/
   filename _infile url 
      "http://research.stlouisfed.org/fred2/data/UNRATE.txt" 
      debug lrecl=100;

   data raw;
      infile _infile missover firstobs = 22;
      format date date9.;
      input @1 date yymmdd10. @13 value 4.1;
   run;

   data _null_;
      set raw end = eof;
      if eof then do;
         interval = intck('month', date, input("&predtime", monyy7.));
         call symput('interval', interval);
         call symput('eodate', date);
         call symput('insert', 'Lastest data:' || put(value, 4.1) || 
                   '% on ' || put(date, monyy7.));
      end;
   run;

   %if %eval(&interval) le 0 %then %do; 
      %put ERROR: Predicted time must be greater than latest time FED posts data; 
      %goto finish; 
   %end;

   ods select none;
   proc forecast data = raw out = _predbyfc outfull
      method = stepar lead = &interval interval = month;
      id date;
      var value;
   run;

   proc arima data = raw;
      identify var = value;
      estimate p = 1 q = 12;
      forecast lead = &interval interval = month 
              id = date out = _predbyar;
   quit;
   ods select all; 

   proc sql;
      create table predicted0 as
      select a.date, a.value label = 'Real unemployment rate', 
            a.forecast as predbyar label = 'ARIMA model', 
            b.value as predbyfc label = 'STEPAR model'
      from _predbyar as a,
           _predbyfc (where = (lowcase(_type_) = 'forecast')) as b
      where a.date = b.date
   ;quit;

   data predicted1;
      set predicted0 end = eof;
      if date lt &eodate then call missing(predbyar, predbyfc);
      else if date eq &eodate then do;
         predbyar = value; 
         predbyfc = value;
      end;
      if eof then do;
         call symput('arlast', put(predbyar, 4.2));
         call symput('fclast', put(predbyfc, 4.2));
      end;
   run;
    
   ods html style = harvest;
   proc sgplot data = predicted1;
      where date ge '01jan2006'd;
      series x = date y = value;
      series x = date y = predbyar;
      series x = date y = predbyfc;
      refline &arlast / axis = y labelloc = inside
         label = "&arlast" transparency = 1;
      refline &fclast / axis = y labelloc = inside
         label = "&fclast" transparency = 1;
      xaxis grid label = ' ';
      yaxis grid label = 'Unemployment percentage %'
         values = (4 to 11 by 0.2);
      inset "Prediction ends on &predtime" / position = topright border;
      inset "&insert" / position = bottomright;
   run;
   ods html close;

   proc datasets nolist;
      delete _:;
   quit;

   %finish: ; 
%mend;

****************(2) TESTING STEP******************;
%unemrate(predtime = NOV2012);

****************END OF ALL CODING***************************************;

Proc Fcmp(4): Binomial tree vs. Black-Scholes model

  
The very truth is that SAS has limited financial functions. Thanks to SAS Institute, they finally added some option pricing functions in the base module of SAS 9.2, such as Black-Scholes put/call functions, Garman-Kohlhagen put/call functions, etc. Thus, the number of financial functions in the SAS toolbox adds up to more than 20 now.

Functions made easy by Proc Fcmp. In the finance’s brave new world of functions, financial function is ammo in a war. A unique and resourceful stockpile of functions is desired. Previously, in SAS, macro may seem like a reasonable substitute. However, when we deal with situations like evaluating y=h(x1)*f(x2), macro is feeble. The new procedure debuted in SAS 9.2, Proc Fcmp, provides us an easy way to create and accumulate functions, which will benefit people who are struggling to use SAS to solve finance problems in Data Step.

Proc Fcmp is used for a Binomial tree function of European call option pricing in this example. Once set up, it is very handy to use a function in the common Data Step: the function just acts like an inherent one. In the example, the exercise price ‘E’ is 50, the time to maturation 't' is 5 months, the share price ‘S’ is 50 again, the risk-free interest rate ‘r’ during t is 0.05, and the volatility ‘sigma’ is 0.3. Then the SAS native Black-Scholes option call function ‘blkshclprc’ and the newly ‘manufactured’ Binomial tree function ‘Eurocall’ powered by Proc Fcmp are compared. With the increasing of the layers ‘n’ of the Binomial tree model or the expanding of the tree branches, the option price fluctuates and gets closer to the price by Black-Scholes model. Finally it may converge to the ‘correct’ price. The results show that more steps may provide more accurate result for a Binomial tree model.

Everything in Proc Fcmp is encapsulated. Typically in designing a macro, we concern about the local or global variables. While dealing with a function by Proc Fcmp, it is argument in and argument out: no variable in a function would leak. All arrays in Proc Fcmp are also temporary unless they are indicated as output. Those features are friendly in making a workable function rapidly.

Reference: 1. Paolo Brandimarte. Numerical Methods in Finance and Economics: A MATLAB-Based Introduction. Wiley-Interscience, 2006.
2. SAS 9.2 Language Reference: Dictionary, Third Edition. SAS Publishing. 2009.

proc fcmp outlib = myfunc.finance.price;
    function Eurocall(E, t, F, r, sigma, N);
        deltaT=T/N;
        u=exp(sigma*sqrt(deltaT));
        d=1/u;
        p=(exp(r*deltaT)-d)/(u-d);
        array tree[1]/ nosymbols;
        call dynamic_array(tree, n+1, n+1);
        call zeromatrix(tree);
        
        do i=0 to N;
            tree[i+1,N+1]=max(0 , E*(u**i)*(d**(N-i)) - F);
        end;
        do j=(N-1) to 0 by -1;
            do  i=0 to j by 1;
                tree[i+1,j+1] = exp(-r*deltaT)* 
                        (p * tree[i+2,j+2] + (1-p) * tree[i+1,j+2]);
            end;  
        end;
        price = tree[1,1];
        return(price);
    endsub;
run;


*****(2)Use Binomial tree model and Black-Scholes model functions *****;
options cmplib = (myfunc.finance);
data test;
    BSprice=blkshclprc(50, 5/12, 50, 0.05, 0.3);
    do n=1 to 100;
        Treeprice=eurocall(50, 5/12, 50, 0.05, 0.3, n);
        output;
    end;
run;

***********(3)Display the comparision between the two functions***************;
proc sgplot data=test;
    title 'The comparison between Black-Sholes model and Binomial tree model';
    needle x=n y=Treeprice/baseline=4;
    series x=n y=BSprice/ lineattrs=(color=red);
    yaxis label='Option price';
run;
***************END***************TEST PASSED 12DEC2010**************;

Proc Fcmp(2): a subroutine for Binomial-CRR model

Problems: Quote for six-month American style euro currency options on plain vanilla, Max[S-K,0]and 〖Max[S-K,0]〗^0.5. Exchange rate S_0=$1.3721 /euro
Six-month continuously compounded inter-bank rates: r=0.4472%,r_f=1.2840%.
Assumptions:The exchange rate for euro follows an iid log normal price changes and volatility is constant.
Methodology:Binomial Model is used to price American currency options on euros.
We calculate the payoffs at time T and discount payoffs to the prior time step. Under the risk neutral probability measure,
c_(t-1)=(q×c_u+(1-q)×c_d)/R
Since these two options are American styles, we need to check for optimal early exercise at each node of the binomial tree. For these two currency options, we check Max[S-K,c_(t-1),0] and Max[〖Max[S-K,0]〗^0.5,c_(t-1) ] . Matlab is the software used to implement the binomial model.
--Parameters:1. Time steps n
As the number of time steps n increases, we can apply CRR model and the binomial model approaches real world price changes. We choose n=80,h=T/n=0.5/80=0.00625.
2. u and dWe choose CRR model to define u and d for binomial model.
u=e^(σ√h) ; d=1/u=e^(-σ√h)
Where h is the length of the binomial times step and σ is the annualized log price change volatility.
3. Annualized log price change volatility σThe daily log changes and daily squared log changes for two year exchange rates from 11/5/2008 – 11/5/2010 are as follows. We consider the volatility to be constant since 05/01/2009. Thus, we choose the historical prices from 05/01/2009-11/5/2010 and apply the volatility of the daily log changes as an estimate. Then the annualized log price change volatility equals the square root of the trading days in one year (252 days) times the daily log price change volatility.
σ=√252×σ_daily
4. Risk Neutral Measure QOptions on currencies can be regarded as an asset providing a yield at the foreign risk-free rate of interest. Thus, the risk neutral probability measure Q :
q=(e^((r-r_f ) )-d)/(u-d) ;
1-q= 〖u- e〗^((r-r_f ) )/(u-d)
5. Strike Price KSet strike price K from $1.3000/euro to $1.5000/euro with $0.005 per euro increments.
Reference: 1. John C. Hull.Options, Futures and Other Derivatives, 7th edition. Prentice Hall. 2008.
2. Base SAS 9.2 Procedures Guide. SAS Publishing. 2009

**********************AUTHOR(DAPANGMAO)----HCHAO8@GMAIL.COM***********************************;
****************(1) CONSTRUCT € TO $ EXCHANGE RATIO VECTOR******************;
data vector;
attrib value informat=dollar10.4 format=dollar10.4 ;
StrikeS=1.3000;  
StrikeE=1.5000;
Increment=0.005;
do value=StrikeS to StrikeE by  Increment;
drop StrikeS StrikeE Increment;
output;
end;
run;

**************(2) BUILD THE FUNCTION TO EVALUATE OPTION PRICES********;
proc fcmp outlib = myfunc.finance.subrout;
subroutine mysubr(T, n, r, rf, s0, sigma, c1[*], c2[*]);
/*Two vectors are output arguments*/
outargs c1, c2;
/*Inside calculation from input arguments*/
dt=T/n;
length=dim(c1);
u=exp(sigma*sqrt(dt));
d=1/u;
q=(exp((r-rf)*dt)-d)/(u-d);
/*Announce 4 arrays -- 1 vector and 3 matrixes */
array k[1]/ nosymbols;
array s[1] /nosymbols;
array x1[1] /nosymbols;
array x2[1] /nosymbols;
n=n+1;
/*The sizes of the arrays are specified*/
call dynamic_array(s, n, n);
call dynamic_array(x1, n, n);
call dynamic_array(x2, n, n);
call dynamic_array(k, length);
/*Read the exchange ratio into function*/
rc=read_array('vector', k);
/*Assign values to S matrix */
call zeromatrix(s);
S[1,1]=S0;
do _col=2 to n ;
do _row=1 to n;
S[_row,_col]=S0*u**(_col-_row)*d**(_row-1);
if _row gt _col then S[_row, _col]=0;
end;
end;
/*Generate final option vectors */
do i=1 to length;
x=k[i];
call zeromatrix (x1);
call zeromatrix (x2);
do j=1 to n;
x1[j,n]=max(S[j,n]-x,0);   
x2[j,n]=max(S[j,n]-x,0)**0.5;     
end;
do _col=(n-1) to 1 by -1;
do _row=1 to (n-1) by 1;
h=exp(-r*dt)*(q*x1[_row,_col+1]+(1-q)*x1[_row+1,_col+1]);
h2=exp(-r*dt)*(q*x2[_row,_col+1]+(1-q)*x2[_row+1,_col+1]);
x1[_row,_col]=max(S[_row,_col]-X,h,0);    
x2[_row,_col]=max( max(S[_row,_col]-x ,0) **0.5, h2);  
end;           
c1[i]=x1[1,1];
c2[i]=x2[1,1];
end;
end;
endsub;
run;
quit;

***********SOME INITIAL VALUES*****************;
/*T=1/2;*/
/*n=80;*/
/*r=0.004472;    */
/*rf=0.01284;*/
/*S0=1.3721;*/
/*Sigma=0.1074*/

**************(3) MEASURE THE LENGTH OF PRICING VECTOR**************;
proc sql;
select count(*) into: vecnum from vector;
quit;

**************(4) USE THE SUBROUTINE TO GENERATE TWO VECTORS********************;
options cmplib = (myfunc.finance);
data final;
array c1[&vecnum] _temporary_;
array c2[&vecnum] _temporary_;
call mysubr(0.5, 80, 0.004472, 0.01284, 1.3721, 0.1074, c1, c2);   /*subroutine mysubr(T, n, r, rf, s0, sigma, c1[*], c2[*]);*/
do i=1 to dim(c1);
c1value= c1[i];
c2value=c2[i];
output;
end;
run;
****************TEST PASSED ----------- END****************************;