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