c*********************************************************************** c Finds maximum likelihood estimates of parameters of symmetric, c mean-reverting, correlated diffusions from discrete time c observations. State variables assumed to follow processes dx_i(t): c c dx_i = k ( u - x_i ) dt + s ( p^1/2 dz_0 + (1-p)^1/2 dz_i ) c c where k, u are mean reversion parameters, s is common total volatilty c and p is common pairwise instantaneous correlation coefficient. c Processes z_0 and z_i, i=1..n, are independent standard Weiner c processes. c State variables are assumed observed at IT discrete intervals of c length H. c c Inputs: N(IT) vector no. of correlated diffusions each t c IT no. of observation dates c H length of time between observations c X(*) state variable levels at start of each interval c concatenated: t=1 vector, t=2 vector, etc. c XH(*) state variable levels time H later of same c IFLAG integer flag selecting options: set to c IMr + 2 * ICov + 4 * IGenZ + 8 * INsame c where c IMr = 0 (no drift) or 1 (mean-reverting) c ICov = 0 (skip COV) or 1 (return COV) c IGenZ = 0 (skip ZHAT) or 1 (return ZHAT) c INsame = 0 (N's vary) or 1 (N(1) applies to all) c c Outputs: K estimated mean reversion coefficient c U estimated mean reversion target c SIG estimated instantaneous volatility c RHO estimated common correlation coeffient c LMAX maximized 2 x log-likelihood value c ZHAT() estimated common factor series c COV lower triangular part of est. cov. matrix c c Parameters: MAXT max number of observation dates (change to suit) c c Otherwise no limitations on size of data arrays c c Author: R. Jones 4 June 1999 c add estim common factor series as output 24/8/99 c add estim cov matrix of parameters 24/9/99 c change to BHH covariance matrix estimate 9/10/99 c change to permit incomplete data series 13/10/99 c*********************************************************************** subroutine & MAXLIKE( N, IT, H, X, XH, IFLAG, K, U, SIG, RHO, LMAX, ZHAT, COV) implicit double precision (A-H, K-L, O-Z) logical MREVERT, GENCOV, GENZHAT, NSAME external LAMBDA parameter (MAXT = 500) ! max no. of time points parameter (EPS = 1d-2) dimension X(*), XH(*), ZHAT(IT), COV(10), N(IT) dimension xe(MAXT), xhe(MAXT), xx(MAXT), xxh(MAXT), & xhxh(MAXT), g(MAXT), NN(MAXT) common / LLCOM / XFLAG, DAT(18), xe, xhe, xx, xxh, xhxh, g, NN equivalence ( DAT( 1), XT ), & ( DAT( 2), XN ), & ( DAT( 3), Sxe ), & ( DAT( 4), Sxx ), & ( DAT( 5), Sxhe ), & ( DAT( 6), Sxhxh ), & ( DAT( 7), Sxxh ), & ( DAT( 8), Skxe ), & ( DAT( 9), Skxhe ), & ( DAT(10), Skxeex ), & ( DAT(11), Skxeexh ), & ( DAT(12), Skxheexh), & ( DAT(13), Skn ) if ( IT .gt. MAXT ) then print*,'MAXLIKE: Observation dates exceeds MAXT' return endif c Determine which options apply if ( IFLAG .gt. 15 .or. IFLAG .lt. 0 ) then print*, 'MAXLIKE.FOR: Inadmissible Iflag = ',IFLAG return endif II = IFLAG if ( mod(II,2) .eq. 1 ) then MREVERT = .true. else MREVERT = .false. endif II = II / 2 if ( mod(II,2) .eq. 1 ) then GENCOV = .true. else GENCOV = .false. endif II = II / 2 if ( mod(II,2) .eq. 1 ) then GENZHAT = .true. else GENZHAT = .false. endif II = II / 2 if ( mod(II,2) .eq. 1 ) then NSAME = .true. else NSAME = .false. endif c Copy some stuff to COMMON XT = dble(IT) if ( NSAME ) then do I = 1, IT NN(I) = N(1) enddo else do I = 1, IT NN(I) = N(I) enddo endif c Accumulate sample moments of the data do I = 2, 13 DAT(I) = 0d0 enddo c For each observation date IX = 0 ! position in X() vector XN = 0 ! total no. of observations do J = 1, IT IN = NN(J) xx(J) = 0d0 xe(J) = 0d0 xxh(J) = 0d0 xhe(J) = 0d0 xhxh(J) = 0d0 XN = XN + IN do I = 1, IN y = X (IX+I) yh = XH(IX+I) xx(J) = xx(J) + y * y xxh(J) = xxh(J) + y * yh xhxh(J) = xhxh(J) + yh * yh xe(J) = xe(J) + y xhe(J) = xhe(J) + yh enddo IX = IX + IN enddo c Accumulate total moments not depending on RHO Sxx = 0d0 Sxxh = 0d0 Sxhxh = 0d0 Sxe = 0d0 Sxhe = 0d0 do J = 1, IT Sxx = Sxx + xx(J) Sxxh = Sxxh + xxh(J) Sxhxh = Sxhxh + xhxh(J) Sxe = Sxe + xe(J) Sxhe = Sxhe + xhe(J) enddo c Determine case if ( MREVERT ) then ! mean-reverting XFLAG = 1d0 else ! zero drift A = 1d0 B = 0d0 S = (Sxhxh - 2d0 * Sxxh + Sxx) / XN XFLAG = 0d0 endif c Find rho maximizing concentrated likelihood function ax = .00d0 ! initial guess range bx = +.03d0 call MNBRAK(ax,bx,cx,fa,fb,fc,LAMBDA) ! bracket maximum tol = 1d-8 LMAX = - BRENT(ax,bx,cx,LAMBDA,tol,RHOD) ! find maximum c Compute remaining parameters a,b,s of discrete-time model: if ( MREVERT ) then RHO = RHOD TOP = Skn * Sxxh - RHO * Skn * Skxeexh - (1d0-RHO)*Skxe*Skxhe BOT = Skn * Sxx - RHO * Skn * Skxeex - (1d0-RHO)*Skxe**2 A = TOP / BOT B = (Skxhe - A * Skxe) / Skn TOP = Sxhxh + A*A*Sxx + B*B*XN - 2d0*(A*Sxxh + B*Sxhe -A*B*Sxe) S = TOP / XN endif c Convert to continuous-time model parameters: if ( .not. MREVERT ) then K = 0d0 U = 0d0 SIG = sqrt( S / H ) elseif ( MREVERT ) then if ( A .gt. 0d0 ) then K = - log(A) / H else K = 1d1 / H A = exp( - K * H ) print*,' Maxlike: warning A < 0 ... kappa = 10 output' endif U = B / ( 1d0 - A ) SIG = sqrt( 2d0 * S * log(A) / ( H * (A**2 - 1d0 ) ) ) endif RHO = RHOD c Get true 2 * log-likelihood function value: LMAX = LMAX - XN * log(4d0 * datan(1d0)) c Generate estimated common factor movements: if ( GENZHAT ) then if ( S * RHO .lt. 1d-10 ) then ! indeterminate or RHO < 0 do J = 1, IT ZHAT(J) = 0d0 enddo Sze = 0d0 Sze2 = 0d0 else ! positive RHO Sze = 0d0 Sze2 = 0d0 do J = 1, IT IN = NN(J) scale = 1d0 / ( IN * sqrt(S * RHO / H ) ) ze = ( xhe(J) - A * xe(J) - IN * B ) * scale ZHAT(J) = ze Sze = Sze + ze Sze2 = Sze2 + ze * ze enddo endif endif c Generate parameter cov estimate using Berndt/Hall/Hausman: if ( GENCOV ) then do I = 1, 10 COV(I) = 0d0 enddo do J = 1, IT c Generate numerical first partials of log Likelihood D = EPS * SIG + 1d-5 DLDS = (LOGL(K,U,SIG+D,RHO,H,J) & - LOGL(K,U,SIG-D,RHO,H,J)) / (2d0 * D) D = EPS * RHO + 1d-5 DLDR = (LOGL(K,U,SIG,RHO+D,H,J) & - LOGL(K,U,SIG,RHO-D,H,J)) / (2d0 * D) if ( MREVERT ) then D = EPS * K + 1d-5 DLDK = (LOGL(K+D,U,SIG,RHO,H,J) & - LOGL(K-D,U,SIG,RHO,H,J)) / (2d0 * D) D = EPS * U + 1d-5 DLDU = (LOGL(K,U+D,SIG,RHO,H,J) & - LOGL(K,U-D,SIG,RHO,H,J)) / (2d0 * D) endif c Accumulate cross product of partials into COV order S,R,K,U COV(1) = COV(1) + DLDS * DLDS COV(2) = COV(2) + DLDS * DLDR COV(3) = COV(3) + DLDR * DLDR if ( MREVERT ) then COV(4) = COV(4) + DLDS * DLDK COV(5) = COV(5) + DLDR * DLDK COV(6) = COV(6) + DLDK * DLDK COV(7) = COV(7) + DLDS * DLDU COV(8) = COV(8) + DLDR * DLDU COV(9) = COV(9) + DLDK * DLDU COV(10)= COV(10) + DLDU * DLDU endif enddo c Invert information matrix estimate to get COV matrix if ( .not. MREVERT ) then call RSINV(COV,2,IER) elseif ( MREVERT ) then call RSINV(COV,4,IER) endif if (IER .eq. -1) then ! return 0 if not + def do i = 1, 10 COV(i) = 0d0 enddo endif endif ! end of GENCOV block return end c*********************************************************************** c Function calculates SINGLE DAY log-likelihood function for correlated c diffusions with mean reversion from data moments and parameters. c c Arguments: KAP, U, SIG, R values of model parameters c H time interval between observations c I observation number c c Common : DAT a vector containing the data moments as follows: c 1 N number of jointly correlated diffusions c 2 T number of discrete time observations c 3 Sum_t x'e c 4 Sum_t x'x c 5 Sum_t xh'e c 6 Sum_t xh'xh c 7 Sum_t x'ee'x c 8 Sum_t xh'ee'xh c 9 Sum_t x'xh c 10 Sum_t x'ee'xh c Author: R. Jones 15/10/99 c*********************************************************************** double precision function LOGL( KAP, U, SIG, R, H, I ) implicit double precision (A-Z) integer NN, MAXT, I parameter (MAXT = 500) dimension xe(MAXT), xhe(MAXT), xx(MAXT), xxh(MAXT), & xhxh(MAXT), g(MAXT), NN(MAXT) common / LLCOM / XFLAG, DAT(18), xe, xhe, xx, xxh, xhxh, g, NN c Transform to discrete time model parameters A = exp( - H * KAP ) B = ( 1d0 - A ) * U if ( KAP .eq. 0d0 ) then S = SIG ** 2 * H else S = SIG ** 2 * ( 1d0 - A*A ) / ( 2d0 * KAP ) endif R1 = 1d0 - R N = dble( NN(I) ) k = 1d0 / ( 1d0 + ( N - 1d0 ) * R ) PI2 = 4d0 * datan(1d0) Sxx = xx(I) Sxe = xe(I) Sxhe = xhe(I) Sxxh = xxh(I) Sxhxh = xhxh(I) L = - N * log(PI2) - N * log(S) - (N-1d0) * log(R1) + log(k) JNK = Sxhxh + A*A * Sxx - 2d0 * A * Sxxh + B*B * R1 * k * N & - R * k * ( Sxhe - A * Sxe ) ** 2 & - 2d0 * B * R1 * k * ( Sxhe - A * Sxe ) L = L - JNK / (R1 * S) LOGL = L / 2d0 return end c*********************************************************************** c Subroutine calculating CONCENTRATED likelihood function for c estimation of symmetric correlated diffusions with mean reversion. c I.e., 2 log Likelihood maximized with respect to c mean reversion parameters a, b, s, and is a function just of rho c and the moments of the data. c c Arguments: RHO a value for the common correlation coefficient c DATX a vector containing the data moments as follows: c 1 N number of jointly correlated diffusions c 2 T number of discrete time observations c 3 Sum_t x'e c 4 Sum_t x'x c 5 Sum_t xh'e c 6 Sum_t xh'xh c 7 Sum_t x'ee'x c 8 Sum_t xh'ee'xh c 9 Sum_t x'xh c 10 Sum_t x'ee'xh c c Notes : x(t) = column n-vector of state values at time t c xh(t) = column n-vector of state values at time t+h c e = column n-vector of 1's c *** -log L is returned so Numerical Recipes minimization c routines can be used c c Author : 14 October 1999 R. Jones c*********************************************************************** double precision function LAMBDA( RHO ) implicit double precision (A-Z) integer IT, NN, MAXT parameter (MAXT = 500) dimension xe(MAXT), xhe(MAXT), xx(MAXT), xxh(MAXT), & xhxh(MAXT), k(MAXT), NN(MAXT) common / LLCOM / XFLAG, DAT(18), xe, xhe, xx, xxh, xhxh, k, NN equivalence ( DAT( 1), T ), & ( DAT( 2), N ), & ( DAT( 3), Sxe ), & ( DAT( 4), Sxx ), & ( DAT( 5), Sxhe ), & ( DAT( 6), Sxhxh ), & ( DAT( 7), Sxxh ), & ( DAT( 8), Skxe ), & ( DAT( 9), Skxhe ), & ( DAT(10), Skxeex ), & ( DAT(11), Skxeexh ), & ( DAT(12), Skxheexh), & ( DAT(13), Skn ) c (Check next condition for applicability with varying NN(t)) c Constrain RHO to be in range [-1,1] by returning bad value: if ( RHO .lt. (1d-10 - 1d0/(NN(1)-1)) & .or. RHO .gt. (1d0 - 1d-10) ) then LAMBDA = 1d20 return endif Slogz = 0d0 do IT = 1, nint(T) ! needed for both g = 1d0 + ( NN(IT) - 1 ) * RHO Slogz = Slogz + log(g) k(IT) = 1d0 / g enddo if ( XFLAG .eq. 0d0 ) then ! no drift case Skyeey = 0d0 do IT = 1, nint(T) Skyeey = Skyeey + k(IT) * (xhe(IT) - xe(IT))**2 enddo S = (Sxhxh - 2d0 * Sxxh + Sxx) / N L = (N-T) * log(1d0 - RHO) + N * log(S) + Slogz & + N / (1d0 - RHO) - RHO * Skyeey / ((1d0 - RHO) * S) LAMBDA = L return elseif ( XFLAG .eq. 1d0 ) then ! mean reverting Skn = 0d0 Skxe = 0d0 Skxhe = 0d0 Skxeexh = 0d0 Skxeex = 0d0 Skxheexh = 0d0 do IT = 1, nint(T) kt = k(IT) xet = xe(IT) xhet = xhe(IT) Skn = Skn + kt * NN(IT) Skxe = Skxe + kt * xet Skxhe = Skxhe + kt * xhet Skxeexh = Skxeexh + kt * xet * xhet Skxeex = Skxeex + kt * xet ** 2 Skxheexh = Skxheexh + kt * xhet ** 2 enddo TOP = Skn * Sxxh - RHO * Skn * Skxeexh - (1d0-RHO)*Skxe*Skxhe BOT = Skn * Sxx - RHO * Skn * Skxeex - (1d0-RHO)*Skxe**2 A = TOP / BOT B = (Skxhe - A * Skxe) / Skn TOP = Sxhxh + A*A*Sxx + B*B*N - 2d0*(A*Sxxh + B*Sxhe -A*B*Sxe) S = TOP / N JNK = Sxhxh + A*A*Sxx - 2d0*A*Sxxh + B*B*(1d0-RHO)*Skn & - RHO * (Skxheexh - 2d0 * A * Skxeexh + A*A * Skxeex) & - 2d0 * B * (1d0-RHO) * (Skxhe - A * Skxe) L = (N-T) * log(1d0 - RHO) + N * log(S) + Slogz & + JNK / (1d0-RHO) / S LAMBDA = L return else ! flag set wrong LAMBDA = 0d0 return endif end c*********************************************************************** SUBROUTINE mnbrak(ax,bx,cx,fa,fb,fc,func) DOUBLE PRECISION ax,bx,cx,fa,fb,fc,func,GOLD,GLIMIT,TINY EXTERNAL func PARAMETER (GOLD=1.618034d0, GLIMIT=100.d0, TINY=1.d-20) DOUBLE PRECISION dum,fu,q,r,u,ulim fa=func(ax) fb=func(bx) if(fb.gt.fa)then dum=ax ax=bx bx=dum dum=fb fb=fa fa=dum endif cx=bx+GOLD*(bx-ax) fc=func(cx) 1 if(fb.ge.fc)then r=(bx-ax)*(fb-fc) q=(bx-cx)*(fb-fa) u=bx-((bx-cx)*q-(bx-ax)*r)/(2.d0*sign(max(abs(q-r),TINY),q-r)) ulim=bx+GLIMIT*(cx-bx) if((bx-u)*(u-cx).gt.0.d0)then fu=func(u) if(fu.lt.fc)then ax=bx fa=fb bx=u fb=fu return else if(fu.gt.fb)then cx=u fc=fu return endif u=cx+GOLD*(cx-bx) fu=func(u) else if((cx-u)*(u-ulim).gt.0.d0)then fu=func(u) if(fu.lt.fc)then bx=cx cx=u u=cx+GOLD*(cx-bx) fb=fc fc=fu fu=func(u) endif else if((u-ulim)*(ulim-cx).ge.0.d0)then u=ulim fu=func(u) else u=cx+GOLD*(cx-bx) fu=func(u) endif ax=bx bx=cx cx=u fa=fb fb=fc fc=fu goto 1 endif return END C (C) Copr. 1986-92 Numerical Recipes Software *1`3.d0 c*********************************************************************** FUNCTION brent(ax,bx,cx,f,tol,xmin) INTEGER ITMAX DOUBLE PRECISION brent,ax,bx,cx,tol,xmin,f,CGOLD,ZEPS EXTERNAL f PARAMETER (ITMAX=100,CGOLD=.3819660d0,ZEPS=1.d-18) INTEGER iter DOUBLE PRECISION a,b,d,e,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x *,xm a=min(ax,cx) b=max(ax,cx) v=bx w=v x=v e=0.d0 fx=f(x) fv=fx fw=fx do 11 iter=1,ITMAX xm=0.5d0*(a+b) tol1=tol*abs(x)+ZEPS tol2=2.d0*tol1 if(abs(x-xm).le.(tol2-.5d0*(b-a))) goto 3 if(abs(e).gt.tol1) then r=(x-w)*(fx-fv) q=(x-v)*(fx-fw) p=(x-v)*q-(x-w)*r q=2.d0*(q-r) if(q.gt.0.d0) p=-p q=abs(q) etemp=e e=d if(abs(p).ge.abs(.5d0*q*etemp).or.p.le.q*(a-x).or.p.ge.q*(b-x) *) *goto 1 d=p/q u=x+d if(u-a.lt.tol2 .or. b-u.lt.tol2) d=sign(tol1,xm-x) goto 2 endif 1 if(x.ge.xm) then e=a-x else e=b-x endif d=CGOLD*e 2 if(abs(d).ge.tol1) then u=x+d else u=x+sign(tol1,d) endif fu=f(u) if(fu.le.fx) then if(u.ge.x) then a=x else b=x endif v=w fv=fw w=x fw=fx x=u fx=fu else if(u.lt.x) then a=u else b=u endif if(fu.le.fw .or. w.eq.x) then v=w fv=fw w=u fw=fu else if(fu.le.fv .or. v.eq.x .or. v.eq.w) then v=u fv=fu endif endif 11 continue pause 'brent exceed maximum iterations' 3 xmin=x brent=fx return END C (C) Copr. 1986-92 Numerical Recipes Software *1`3.d0 C**** SUBROUTINE TO INVERT REAL SYMMETRIC MATRIX: IBM SUBROUTINE LIBRARY SUBROUTINE RSINV(AMAT,N,IER) IMPLICIT REAL*8(A-H,O-Z) DIMENSION AMAT(1) CALL RMFSD(AMAT,N,IER) IF(IER) 9,1,1 1 IPIV = N*(N+1)/2 IND = IPIV DO 6 I = 1,N DIN = 1.D0/AMAT(IPIV) AMAT(IPIV) = DIN MIN = N KEND = I - 1 LANF = N - KEND IF(KEND) 5,5,2 2 J = IND DO 4 K = 1,KEND WORK = 0.D0 MIN = MIN - 1 LHOR = IPIV LVER = J DO 3 L = LANF,MIN LVER = LVER + 1 LHOR = LHOR + L 3 WORK = WORK + AMAT(LVER)*AMAT(LHOR) AMAT(J) = -WORK*DIN 4 J = J - MIN 5 IPIV = IPIV - MIN 6 IND = IND - 1 DO 8 I = 1,N IPIV = IPIV + I J = IPIV DO 8 K = I,N WORK = 0.D0 LHOR = J DO 7 L = K,N LVER = LHOR + K - I WORK = WORK + AMAT(LHOR)*AMAT(LVER) 7 LHOR = LHOR + L AMAT(J) = WORK 8 J = J + K 9 RETURN END SUBROUTINE RMFSD(AMAT,N,IER) IMPLICIT REAL*8(A-D,F-H,O-S,U-Z) DIMENSION AMAT(1) EPS = 0.0 IF(N-1) 12,1,1 1 IER = 0 KPIV = 0 DO 11 K = 1,N KPIV = KPIV + K IND = KPIV LEND = K - 1 TOL = ABS(EPS*SNGL(AMAT(KPIV))) DO 11 I = K,N DSUM = 0.D0 IF(LEND) 2,4,2 2 DO 3 L = 1,LEND LANF = KPIV - L LIND = IND - L 3 DSUM = DSUM +AMAT(LANF)*AMAT(LIND) 4 DSUM = AMAT(IND) - DSUM IF(I-K) 10,5,10 5 IF(SNGL(DSUM)-TOL) 6,6,9 6 IF(DSUM) 12,12,7 7 IF(IER) 8,8,9 8 IER = K - 1 9 DPIV = DSQRT(DSUM) AMAT(KPIV) = DPIV DPIV = 1.D0/DPIV GO TO 11 10 AMAT(IND) =DSUM*DPIV 11 IND = IND+ I RETURN 12 IER = -1 RETURN END