Translate GAUSS code to R code

Joined
10/12/11
Messages
28
Points
11
ATTENTION: THIS IS NOT VBA. IT'S GAUSS-CODE. Sry for caps.

Hi Folks,

I have some VBA code, that I wish to "translate" into R code. I understand R, but is useless in VBA. Anybody care to help? Thanks.

Code:
/* CCRgDer.txt Program to perform derivative calculations for the risk parameters
for a Contingent Capital Bond where the conversion price is the trigger price.
Derivatives are computed for different values of initial capital to deposits, x.
It uses the same sequence of random numbers to carry out this exercise. */

ns = 100000; /* Number of simulations */
dt = 1/250; /* Length of each subperiod (e.g., trading day) */

/* Bank Assumptions */
x0 = seqa(1.065,0.005,18); /* Starting ratio of assets to deposits */

np = rows(x0);
/* sga = 0.02; */ /* Brownian motion volatility of bank assets */
rar = -0.2; /* Correlation between asset and interest rate Brownian motions */
/* lm = 1; */ /* Annual frequency of jumps */
/* muy = -0.01; */ /* mean log jump size */
/* sgy = 0.02; */ /* Standard deviation of jump size */
g = 0.5; /* Deposit growth mean reversion */
xs = 1.10; /* Target asset to deposit ratio */

/* Contingent Capital Assumptions */
b0 = 0.04; /* Starting ratio of contingent capital to deposits */
p = 0.9; /* Conversion premium or discount */
xcc = 1.02 + b0*(1-p); /* (Asset-p*bt)/deposit or converted equity/deposit level at which contingent capital converts */
tc = 5; /* Maturity in years of contingent capital */
ff = 1; /* 0; */ /* Dummy that equals 0 if coupon is fixed, equals 1 if the coupon is floating */
c = 0.032675; /* Yield or spread for a starting value of x=1.1 and parameters equal to benchmark */

/* CIR Term Structure Assumptions */
r0 = 0.035; /* Beginning real interest rate */
kr = 0.114; /* Risk-neutral mean reversion */
rb = 0.069; /* Risk-neutral steady state */
sgr = 0.0700; /* Standard deviation of short-rate changes */

mxp = tc/dt; /* Maximum number of periods for bonds */

wr = sqrt(1-rar^2); /* Weight on idionsycratic risk for interest rate shock */
/* k = exp(muy+0.5*sgy^2) - 1; */ /* Expected jump size */

/* Calculate risk-free CIR coupon bond yield */
i = 1;
dsum = 0;
thet = sqrt(kr^2+2*sgr^2);
do until i >mxp;
  t=i*dt;
  ati = ((2*thet*exp((thet+kr)*t/2))/((thet+kr)*(exp(thet*t)-1)+2*thet))^(2*rb*kr/sgr^2);
  bti = 2*(exp(thet*t)-1)/((thet+kr)*(exp(thet*t)-1)+2*thet);
  dsum = dsum + ati*exp(-bti*r0)*dt;
  i = i + 1;
endo;
atc = ((2*thet*exp((thet+kr)*tc/2))/((thet+kr)*(exp(thet*tc)-1)+2*thet))^(2*rb*kr/sgr^2);
btc = 2*(exp(thet*tc)-1)/((thet+kr)*(exp(thet*tc)-1)+2*thet);
ccir = (1 - atc*exp(-btc*r0))/dsum; /* 0.05; */

output file = c:\doc\rbc\CCTgRkDrFlb04x02p9.out reset;
outwidth 240;

nps = 5; /* Number of different parameter sets */

vs = zeros(np,nps); /* Average value of risk-neutral coupon bond payoff (Each iteration do vs = vs + val/ns) */

" ";
" x0 = ";;x0;
" ";
" ";
" Risk-free CIR yield ";;ccir;

i = 1;
do until i > ns; /* Do loop for each simulation */

  if fmod(i,1000) == 0;
    " Simulation = ";;i;
  endif;

  ep = rndn(mxp,3); /* Three normally-distributed random numbers for each potential period. */
  unif = rndu(mxp,1); /* .<(lm*dt); */ /* A Poisson count for each period having probability lm*dt */

  ic = 1; /* Do loop for a given coupon and simulation */
  do until ic > nps;

    if ic < 1.5;
      lm = 1;
      psn = unif.<(lm*dt);
      sga = 0.02;
      muy = -0.01;
      sgy = 0.02;
      k = exp(muy+0.5*sgy^2) - 1;

    elseif ic < 2.5;
      lm = 1.25;
      psn = unif.<(lm*dt);
      sga = 0.02;
      muy = -0.01;
      sgy = 0.02;
      k = exp(muy+0.5*sgy^2) - 1;

    elseif ic < 3.5;
      lm = 1;
      psn = unif.<(lm*dt);
      sga = 0.025;
      muy = -0.01;
      sgy = 0.02;
      k = exp(muy+0.5*sgy^2) - 1;

    elseif ic < 4.5;
      lm = 1;
      psn = unif.<(lm*dt);
      sga = 0.020;
      muy = -0.0125;
      sgy = 0.02;
      k = exp(muy+0.5*sgy^2) - 1;

    else;
      lm = 1;
      psn = unif.<(lm*dt);
      sga = 0.020;
      muy = -0.01;
      sgy = 0.025;
      k = exp(muy+0.5*sgy^2) - 1;

    endif;

    ix = 1; /* Do loop for a given x value, parameter set, and simulation */
    do until ix > np;
      x = x0[ix,1]; /* Initialize asset to deposit ratio */

      t = 0;
      j=1; /* j is index for time, t */
      val = 0;
      r = r0; /* Initialize interest rate */
      bt = b0; /* Initialize contingent capital par to deposit ratio */
      dsf = 1; /* Initialize discount factor */

      do until t > (tc-dt+0.000001); /* Do loop for each period over length of bond maturity */

        /* Compute deposit credit spread */
        d1 = (ln(x)+muy)/sgy;
        d2 = d1+sgy;
        h = lm*(cdfn(-d1)-x*exp(muy+0.5*sgy^2)*cdfn(-d2));
        if h < 0;
          h = 0;
        endif;
        /* " t = ";;t;;" x = ";;x;;" h = ";;h; */
        /* Update x for next period */
        x = x*exp((r-lm*k - (r+h+(c+ff*r)*bt)/x - g*(x-xs) - 0.5*sga^2)*dt + sga*sqrt(dt)*ep[j,1] + psn[j,1]*(muy + sgy*ep[j,3]));
        dsf = dsf*exp(-r*dt);

        if x > (xcc + p*bt); /* Continuation state */
          val = val + dsf*(c+ff*r)*b0*dt; /* Discounted end of period coupon added to value. */

          /* if t > maxt;
            maxt = t;
            " maxt = ";;maxt;
          endif; */

          if t > (tc - dt - 0.000001);
            val = val + dsf*b0; /* Payment of final principal */
            /* " Maturity Reached"; */
          endif;

          r = rb*kr*dt + r*(1-kr*dt) + sgr*sqrt(r*dt)*(rar*ep[j,1]+wr*ep[j,2]); /* Update Interest Rate for next period */
          bt = bt*exp(-g*(x-xs)*dt); /* Update cc par to deposit ratio */

          t = t + dt;
          j = j + 1;

        elseif x > 1; /* Conversion without Failure */
          val = val + dsf*(x-1)*p*b0/(xcc+p*bt-1); /* Check if p < 1 case makes sense in terms of number of shares issued. */
          t = tc + 1;

        else; /* Failure: Value not incremented */
          t = tc + 1;

        endif;

      endo; /* Endo for time period */

      vs[ix,ic] = vs[ix,ic] + val/ns;

      ix = ix + 1;
    endo; /* Endo for x value */

    ic = ic + 1;
  endo; /* Endo for parameter set */

  /* " Simulation = ";;i; */
  i = i + 1;

endo; /* Endo for each simulation */

" ";
" vs = ";
vs;
" ";
" ";
pdif = vs/b0;
pdif = ln(pdif);
" pdif = ";;
pdif;

output off;
 
That reinforces the correctness of my idea to invent something like Google Translate for programming languages.
 
Haha.. thought about this like 1E99999 times!
OMG.. the time one would safe... not to mention the "lost-in-translation"-issues would probably expand by a gazillion :)
 
But srsly - if anyone know about programmers that can me hired by the hour for jobs like these, I'd be willing to wire some $$ for a R-program (that works). Right now this is an obstacle, I don't have time to deal with. Please let me know if anyone know anyone who might be interested in this fairly easy task. Thanks.
 
Just got a note from cgorac informing me that this is Gauss, not VBA, code. (Thanks). I'm not that familiar with VBA and not familiar with Gauss at all, so I made my assumption based on the graphs associated with the program (see attached).

I will do a "re-post" shortly, and describe what I need. Sorry for this inconvenience.

Penn4.png
 
If that code was more organized with indentations to line up the ifs and the do (for) loops, I'd be able to give you a program in short order.
 
Typing gg=(shift+g) in VIM will autoindent it, assuming your VIM has an indent package for the language installed.
 
Just downloaded MacVim. I have absolutely no idea of how to operate it. Found a page stating that I should use "set ai"... but how do you execute commands?
 
If that code was more organized with indentations to line up the ifs and the do (for) loops, I'd be able to give you a program in short order.
Okay, we'll just stay in this thread. No need to "spam" the forum more than necessary. Got an "indented"-version from cgorac. See attached.
I updated the first post for OP with the indented version. Given that GeSHI does not support GAUSS, I wrap them as Matlab code. Now, let's see IlyaKEightSix does his R thing.
 
Andy Nguyen, thanks!

To change things a little, I will step out of my comfort zone and show you some of my "ugly-a**-R-code".

This part yields Pennacchi's figure 2 (also attached).

If anyone has some easy-to-implement optimization, please let me know. Right now my program takes 15+ hours before it yield some "pretty" (interpretable) graphs.

And let me just elaborate on what it is, I need: So, basically first I need to expand my program so that it graphs the value of equity plotted against capital-to-deposit. From there on, I need several editions of that particular program. One that plot the changes in value of equity when the jump-intensity is adjusted, and one that plot the changes in value of equity when the average jump-size is adjusted etc..

Thanks to everyone contributing to this thread.
 

Attachments

Alright folks - I actually manage to do the coding myself. I didn't translate the posted Gauss-code. I read the article, and made my own program. Again, I just want to express my gratitude to everyone who have spent/"wasted" some thoughts.

However, trees don't grow to the sky. Right know my program takes approx. 24h+ when compiling :-/ Any experts on how to get R to deal "faster" with loads of for-loops?

Thanks,
J
 
For loops suck in R because it has to re-interpret every line inside the loop on every evaluation. Try to vectorize the code as much as possible.

And if that isn't good enough, consider using Revolution-R as it's built against Intel MKL linear algebra libraries which offer about a 10x or more increase in performance over the base BLAS.
 
Back
Top Bottom