C.hr saints.f
C....*...1.........2.........3.........4.........5.........6.........7.*
C    
C     saints.f
C
C     Copyright (C) 1999.
C
C     A. Ronald Gallant, P. O. Box 659, Chapel Hill NC 27514-0659, USA
C
C     Permission to use, copy, modify, and distribute this software and 
C     its documentation for any purpose and without fee is hereby 
C     granted, provided that the above copyright notice appear in all
C     copies and that both that copyright notice and this permission 
C     notice appear in supporting documentation.
C     
C     This software is provided "as is" without any expressed or implied 
C     warranty.
C
C....*...1.........2.........3.........4.........5.........6.........7.*
C.hr saints
C@
*....*...1.........2.........3.........4.........5.........6.........7.*
*
*     saints     9/29/99
*
*     purpose
*     in the esaints model, factors follow the SDE
*
*       dy = [b+l0+(A+L1)y]dt + [(Sig)**.5]dW,
*          = [Phi0+(Phi1)y]dt + [(Sig)**.5]dW,
*    
*     and the short rate is determined by 
*
*       r = alpha + y'Ry
*
*     where y and W have dimension d.
*
*     given factors y and parameters b, A, Sig, alpha, and R, this 
*     routine computes the bond yields
*
*        yld(tau) = - {eta(tau) + [gam(tau)]'y + y'[Bgam(tau)]y} / tau
*
*     at maturities tau = h, 2h, 3h, ..., nh, where eta(tau), gam(tau), 
*     and Bgam(tau) are obtained by solving the ODE's 
*     
*       (d/dtau)Bgam(tau) = 2(Bgam)(Sig)(Bgam) + (Bgam)A + A'(Bgam) - R
*       (d/dtau)gam(tau)  = 2(Bgam)(Sig)(gam) + A'(gam) + 2(Bgam)b
*       (d/dtau)eta(tau)  = tr[(Sig)(Bgam)] + (1/2)(gam)'(Sig)(gam)
*                             + (gam)'b - alpha
*
*     with initial condition zero at tau=0. one unit of time is a year;
*     therefore, if h=1.0/12.0 and n=120, then yld(3) is the yield 
*     to maturity of a 3 month bond, yld(12) is the yield to maturity 
*     of a 1 year bond, and yld(120) is the yield to maturity of a 
*     10 year bond.  the derivatives are smooth so h can usually be set
*     to the shortest maturity of interest.  i.e. if the shortest is a
*     3 month bond, then put h=1.0/4.0 and yld(1) contains the yield.
*
*     in the esaints model, bond prices are given by
*
*       P(tau) = exp{ eta(tau) + [gam(tau)]'y + y'[Bgam(tau)]y }
*
*     therefore, P(tau) can be recovered from output using 
*
*       P(tau) = exp { - [tau*yld(tau)] }
*     
*     multiply by 100 to convert yields and the short rate to percent.
*
*     several matrices are patterned, as indicated below, but no attempt
*     is made to conserve storage; they are all stored columnwise. 
*
*     usage
*     call saints(y,d,m,h,n,b,A,Sig,alpha,R,yld)
*
*     libf calls
*     dgmaba, dgmprd
*
*     arguments
*     y      - input factors, a d by m matrix containing factors where
*              y(i,t) is factor i at time t, i=1,...,d and t=1,...,m.
*              input, real*8
*     d      - dimension of y.  also, dimension of b and row and column
*              dimensions of A, Sig, and R.  d must be less than 
*              parameter md which is currently set to 15.
*              input, integer*4
*     h      - time increment for solving the ODE's.
*              input, real*8
*     n      - number increments for which ODE's are to be iterated.   
*              also, the dimension of yld.
*              input, integer*4
*     b      - input vector of dimension d
*              input, real*8
*     A      - input lower triangular negative definite d by d matrix 
*              input, real*8
*     Sig    - input diagonal d by d matrix, diagonal entries positive
*              input, real*8
*     alpha  - non-negative scalar
*              input, real*8
*     R      - symmetric d by d matrix with 1.d0 on the diagonal
*              input, real*8
*     yld    - output yields, an n by m matrix containting computed 
*              bond yields to maturity where yld(tau,t) is maturity 
*              tau at time t, tau=1,...,n and t=1,...,m.
*              output, real*8
*     
*     reference
*     Ahn, Dong-Hyun, Robert F. Dittmar, and A. Ronald Gallant (1999),
*     "The Extended SAINTS Model of the Term Structure of Interest 
*     Rates: Theory and Evidence," Manuscript, Department of Finance,
*     Kenan-Flagler Business School, University of North Carolina, 
*     CB# 3490, Chapel Hill NC 27599-3490.
*
      subroutine saints(y,d,m,h,n,b,A,Sig,alpha,R,yld)
      implicit none

*     arguments

      integer*4 d,m,n

      real*8 y,h,b,A,Sig,alpha,R,yld

      dimension y(d,m),b(d),A(d,d),Sig(d,d),R(d,d),yld(n,m)

*     parameters

      integer*4 md,ml
      integer*4 nout
      parameter (md=15)
      parameter (ml=md*md+md+1)
      parameter (nout=3)

*     local variables

      real*8 h2,h6
      real*8 s0,s1,s2,s3
      real*8 d0,d1,d2,d3
      real*8 w1,w2,w3,w4
      real*8 sum

      dimension s0(ml),s1(ml),s2(ml),s3(ml)
      dimension d0(ml),d1(ml),d2(ml),d3(ml)
      dimension w1(md*md),w2(md*md),w3(md),w4(md)

      integer*4 ibBgam,ibgam,ibeta,iaBgam,iagam,iaeta

      integer*4 i,j,k,l,ij,ji,ii

      intrinsic DFLOAT

*     error traps

      if (d.gt.md) then
        write(nout,'(a,i5)') 'Error, saints, md < d =', d
        stop
      end if

*     compute addresses

      ibBgam = 0
      ibgam  = ibBgam + d*d
      ibeta  = ibgam + d

      iaBgam = ibBgam + 1
      iagam  = ibgam  + 1
      iaeta  = ibeta  + 1
      
      l = iaeta

*     initial conditions

      do i=1,l
        s0(i)=0.d0
      end do

*     fourth order Runge-Kutta iterations

      h2=h/2.d0
      h6=h/6.d0

      do k=1,n

*     step one
*     compute d0, derivative at t, and s1, state at t + h/2
*     s1 = s0 + (h/2)*d0
 
      call dgmaba(s0(iaBgam),Sig,w1,d,d)
      call dgmprd(s0(iaBgam),A,w2,d,d,d)
      do i=1,d
      do j=1,d
       ij = -d + d*j + i
       ji = -d + d*i + j
       d0(ibBgam+ij) = 2.d0*w1(ij) + w2(ij) + w2(ji) - R(i,j)
      end do
      end do
      call dgmprd(s0(iaBgam),Sig,w2,d,d,d)
      call dgmprd(w2,s0(iagam),w1,d,d,1)
      call dgmprd(s0(iagam),A,w3,1,d,d)
      call dgmprd(s0(iaBgam),b,w4,d,d,1)
      sum=-alpha
      do i=1,d
        ii = -d + d*i + i
        sum = sum + w2(ii) + s0(ibgam+i)*b(i)
        d0(ibgam+i) = 2.d0*w1(i) + w3(i) + 2.d0*w4(i)
      end do
      call dgmaba(s0(iagam),Sig,w1,d,1)
      d0(iaeta) = sum + w1(1)/2.d0

      do i=1,l
        s1(i) = s0(i) + h2*d0(i)
      end do

*     step two
*     compute d1, derivative at t + h/2, and s2, state at t + h/2
*     s2 = s0 + (h/2)*d1
 
      call dgmaba(s1(iaBgam),Sig,w1,d,d)
      call dgmprd(s1(iaBgam),A,w2,d,d,d)
      do i=1,d
      do j=1,d
       ij = -d + d*j + i
       ji = -d + d*i + j
       d1(ibBgam+ij) = 2.d0*w1(ij) + w2(ij) + w2(ji) - R(i,j)
      end do
      end do
      call dgmprd(s1(iaBgam),Sig,w2,d,d,d)
      call dgmprd(w2,s1(iagam),w1,d,d,1)
      call dgmprd(s1(iagam),A,w3,1,d,d)
      call dgmprd(s1(iaBgam),b,w4,d,d,1)
      sum=-alpha
      do i=1,d
        ii = -d + d*i + i
        sum = sum + w2(ii) + s1(ibgam+i)*b(i)
        d1(ibgam+i) = 2.d0*w1(i) + w3(i) + 2.d0*w4(i)
      end do
      call dgmaba(s1(iagam),Sig,w1,d,1)
      d1(iaeta) = sum + w1(1)/2.d0

      do i=1,l
        s2(i) = s0(i) + h2*d1(i)
      end do

*     step three
*     compute d2, derivative at t + h/2, and s3, state at t + h
*     s3 = s0 + h*d2
 
      call dgmaba(s2(iaBgam),Sig,w1,d,d)
      call dgmprd(s2(iaBgam),A,w2,d,d,d)
      do i=1,d
      do j=1,d
       ij = -d + d*j + i
       ji = -d + d*i + j
       d2(ibBgam+ij) = 2.d0*w1(ij) + w2(ij) + w2(ji) - R(i,j)
      end do
      end do
      call dgmprd(s2(iaBgam),Sig,w2,d,d,d)
      call dgmprd(w2,s2(iagam),w1,d,d,1)
      call dgmprd(s2(iagam),A,w3,1,d,d)
      call dgmprd(s2(iaBgam),b,w4,d,d,1)
      sum=-alpha
      do i=1,d
        ii = -d + d*i + i
        sum = sum + w2(ii) + s2(ibgam+i)*b(i)
        d2(ibgam+i) = 2.d0*w1(i) + w3(i) + 2.d0*w4(i)
      end do
      call dgmaba(s2(iagam),Sig,w1,d,1)
      d2(iaeta) = sum + w1(1)/2.d0

      do i=1,l
        s3(i) = s0(i) + h*d2(i)
      end do

*     step four
*     compute d3, derivative at t + h, and s4, state at t + h
*     s4 = s0 + h * (d0 + d1 + d2 + d1 + d2 + d3) / 6
*     overwrite s0 by s4
 
      call dgmaba(s3(iaBgam),Sig,w1,d,d)
      call dgmprd(s3(iaBgam),A,w2,d,d,d)
      do i=1,d
      do j=1,d
       ij = -d + d*j + i
       ji = -d + d*i + j
       d3(ibBgam+ij) = 2.d0*w1(ij) + w2(ij) + w2(ji) - R(i,j)
      end do
      end do
      call dgmprd(s3(iaBgam),Sig,w2,d,d,d)
      call dgmprd(w2,s3(iagam),w1,d,d,1)
      call dgmprd(s3(iagam),A,w3,1,d,d)
      call dgmprd(s3(iaBgam),b,w4,d,d,1)
      sum=-alpha
      do i=1,d
        ii = -d + d*i + i
        sum = sum + w2(ii) + s3(ibgam+i)*b(i)
        d3(ibgam+i) = 2.d0*w1(i) + w3(i) + 2.d0*w4(i)
      end do
      call dgmaba(s3(iagam),Sig,w1,d,1)
      d3(iaeta) = sum + w1(1)/2.d0

      do i=1,l
        s0(i) = s0(i) + h6*( d0(i) + d3(i) + 2.d0*( d1(i) + d2(i) ) )
      end do

      do j=1,m

        call dgmaba(y(1,j),s0(iaBgam),w1,d,1)
        sum = w1(1) + s0(iaeta)
        do i=1,d
          sum = sum + y(i,j)*s0(ibgam+i)
        end do
        
        yld(k,j) = - sum / (DFLOAT(k)*h)
  
        end do

      end do
  
      return
      end