/*
** Version 01
** Date : 20/07/2001                                        
** --------------------------------------------------------------------------
** Copyright K.Cuthbertson and D. Nitzsche
** "Financial Engineering:Derivatives and Risk Manangement" - J. Wiley 2001
**
**   The formulae used can be found in Appendix Cht 5.2
**
**   Program to calculate forward rates, term structure and US, T-Bill Futures 
**   prices and graph the term structure of spot and forward rates. 
**
*/

new ; cls ;
format /m1/rd 12,6;  output on ; screen on ;

/*
 " Spot rates, propn p.a.        (r  = r01,r02,r03, ... )               " ;
 " Maturity of spot rates: years (t  = t01,t02,t03, ... )               " ;
 " ONE period forward rates      (f  = f12,f23,f34,.. )                 " ; 
 " Successive time periods       (tij = t12,t23,t34,.. )                 " ; 
 " Forward rates anchored at t=1 (f1 = f12,f13,f14,f15,.. )             " ;
 " Time periods anchored at t=1  (t1 = t12,t13,t14,t15,.. )             " ;
 " rc,fc are continuously compounded rates (propn, p.a.                 " ;
 " T-Bill Futures price (Fprice)                                        " ;
 " Discount rate of futures: simple rate % p.a. (discFut)               " ;
 " IMM index (Qf)                                                       " ; 
 " Spot price of T-Bill: Face value 100 (PTbill)                        " ;
 " Discount rate of spot T-Bill: simple rate % p.a.(discTB)             " ;
 " Price of zero coupon bond paying $1 (pzero = 1/(1+r)^n )             " ;
 */

@ ------------------------------  USER INPUT (See Appendix 5.2) -------------------------------- @
 
r     = { 0.09, 0.10, 0.12 } ;    @ Spot rates, compound p.a.       @
tdays = { 32,   122,   212 } ;    @ Horizon's (days) for spot rates @
year  = 365 ;

@ ----------------------------------- END  USER INPUT --------- -------------------------------- @

t = tdays/year ;

tdifDay = tdays - lag1(tdays) ;
tdifDay = trimr(tdifDay,1,0) ; 
tdifDay = miss(0,0)|tdifDay ;

@ ------------------------------------  Calculate US T-Bill Futures -------------------------------------- @

 {pzero,r,rc,f,fc,fc1,tij,t1,fwdc} = forward(r,t) ;             @ Activate the procedure - see below @
 
  Fprice1   = 100./ (1+f).^(90/365)  ; 
  Fprice2   = 100.*exp( -fc.*90/365) ;
  discFut   = (360/90).*(100 - Fprice1) ;
  Qf        = 100 - discFut ; 

@ ------------------------------------  SPOT T-Bill prices and discount rate -------------------------------- @
    PTbill = 100./(1+r).^t ;
    discTB  =  ( (100 -PTbill)./100  ).* (360/tdays) ;

@ ------------------------------------------------ Print Output  -------------------------------------------- @


?; " Price Zero     spot rate    spot(c.c.)  {r1,f12,f23}  {rc1,fc12,fc23,})  {0,fc12,fc13..} "; pzero~r~rc~f~fc~fc1 ;
?; "    Days        Chg in days  Chg in years                                             "; tdays~tdifDay~tij ;
?; " Forward Pr1   Forward Pr2 Fut.Disc Rate   IMM index                                 "; Fprice1~Fprice2~discFut~Qf ;            
?; " Spot Tbill Pr. Spot rates   Years to Mat.                                            "; PTbill~discTB~t ;



 
@ ---------------------------------------------  MORE USER INPUT HERE ------------------------------------- @
@ ---- Set up a more general term structure and calc. forward rates {f12,f23,f34..} and {f12,f13,f14..} --- @ 
 
t = { 1,2,3,4,5,6,7,8,9,10  } ;         @ Use t in years @

 r = zeros(10,1) ;
 
   r[1,1]  = 0.05 ; 
   r[2,1]  = 0.06 ; 
   r[3,1]  = 0.068 ; 
   r[4,1]  = 0.074 ; 
   r[5,1]  = 0.08 ; 
   r[6,1]  = 0.084 ; 
   r[7,1]  = 0.088 ; 
   r[8,1]  = 0.09 ; 
   r[9,1]  = 0.092 ; 
   r[10,1] = 0.093 ; 


{pzero,r,rc,f,fc,fc1,tij,t1,fwdc} = forward(r,t) ;                        @ Activate the procedure - see below @
 
   f1 = exp(fc1) - 1 ; 

?; " ---------------------------------- Second Output - More Spot and Forward Rates --------------------------- " ;
?;
 " Spot rates(comp. and c.c.) and maturities:ie.(r  = r01,r02,r03, ... ),(t=t01,t02,t03,.. )        ";  r~rc~t ;
 " ONE period forward rates(comp, and c.c): f and fc: (ie. f12,f23,f34,.. )                         ";  f~fc ; 
 " Forward rates( comp and c.c.) and time anchored at t=1 (f1=f12,f13,f14,.. ),(t1=t12,t13,t14,.. ) ";  f1~fc1~t1 ;
?;           
?; " t(yrs)     {r1,fc12,fc23} {0,t12,t13,t14} {0,fc12,fc13,fc14}                                   ";  t~fc~t1~fc1 ;
?;
?; " comparing two methods for calc of {fc12,fc13,fc14,..}                                          "; fc1~fwdc ; 


@ -----------------------------------------  Graph Spot and forward Yield Curve --------------------------------- @ 
 
library pgraph;
graphset;

y = rc~fc1 ;
    title("Spot and forward rates(fc12,fc13,fc14,fc15): Yield Curves") ; 
  xlabel("Maturity, Years");
  ylabel("Spot rates and forward rates: both c.c.");
xy(t,y);


@ ----------------------------------------- END OF PROGRAM ------------------------------------------- @

/* ----------------------------------------------------------------------------------------------------
                    PROCEDURE:  Calculate the forward rates, given the spot rates etc. 
   ---------------------------------------------------------------------------------------------------- */

proc(9) = forward(r,t);
 
 local i, pzero, rc, f, fc, fc1, tij, t1, fwdc, fx  ;


    pzero = 1./ (1 + r).^t;        @ price of zero = discount rates @
    tij  = t - lag1(t) ;
    tij  = trimr(tij,1,0) ;          @ tij is vector of  (r-1) x 1 @


                      
    rc   = ln(1+r) ;               @ spot rate must be contin. comp. for use in forward rate formulae @

    fc   = zeros(rows(r),1);

    i = 1;
     do while i <= rows(r) - 1;
       fc[i+1] =  ( t[i+1].*rc[i+1] - t[i].*rc[i] )./ tij[i] ;
    i = i + 1;
    endo;

    fc[1] = rc[1] ;                           @ fc is now {rc1,fc12,fc23,fc34...} @ 
    f     = exp(fc) - 1 ;                     @ ordinary compound rate @
   
    
    tij = zeros(1,1)|tij ;                  @ tij is now  {0,t12,t23,t34..} ~ (r x 1 vector ) @
    t1 = zeros(rows(r),1) ;               

     i = 1;
      do while i <= rows(r) -1 ;
        t1[i+1] = t1[i] + tij[i+1] ;       @ t1 = cum. sum of tij and gives {0,t12,t13,t14, ..}  @      
     i = i + 1;
     endo;

/*
**  Let  fc1(1) = 0,   fc1(2)= fc12,  fc1(3) = fc13, fc1(4) = fc14 that is the forward rates, ALL starting at t=1. 
**  Let tij(1)=0, tij(2)=t12, tij(3)= t23 etc, 
**  Let t1(1) = 0, t1(2) = t12, t1(3) = t13, t1(4) = t14 etc.
**
**  Then for example, the recursion for fc14 =fc1(4) = [ t1(3)*fc1(3) + tij(4)*fc(4) ] / t1(4)
**  or equivalently t14 x f14 = t13 x f13 + t34 x f34  
*/

    fc1 = zeros(rows(r),1);              
 
    i = 1;
     do while i <= rows(r) - 1;
      fc1[i+1] =  ( t1[i].*fc1[i] + tij[i+1].*fc[i+1] )./ t1[i+1] ;
    i = i + 1;
    endo;                                       @  fc1 = {0, fc12, fc13, fc14 .... } @ 

    fc1[1,1] = miss(0,0) ;                       @ 1st entry is missing value, for graph below @  

/*
** Alternative recursion for calculating the forward rates anchored at t=1 (ie. f12,f13,f14 ..) 
** using only the spot rates.
** The recursion is f12 = (2/1) r2 - (1/1) r1, equiv to 2(r2) = r1 +   f12 
**                  f13 = (3/2) r3 - (1/2) r1  equiv to 3(r3) = r1 + 2(f13) etc.
*/

fwdc = zeros( rows(r),1 ) ; 
fx   = 1 ;                         @ --- fixed anchor for forward rate start point - here=1 year --- @ 

i = 2 ; 
do until i > rows(r) ; 
   fwdc[i,.] = ( t[i,1]./t[i-fx,1] ).*rc[i,1] - ( fx./t[i-fx,1] ).*rc[fx,1] ; 
i = i+1 ; 
endo ;

    retp(pzero,r,rc,f,fc,fc1,tij,t1,fwdc);

endp;

/* ------------------------------------------------------------------------------------------------------------ */
end ;