! 4.4.f90 ! 2024-09-24 ! $Id: 4.4.f90 1.2 2024/12/07 08:34:26 s Exp $ implicit none integer i,j,m,m1,m2,n,n1,n2,nn parameter (m=18,n=9,nn=20, m1=18,n1=9, m2=9,n2=5) double precision face_value double precision ask_price(m1) double precision cpn_rate(m1), cpn(m1) double precision accrued_interest(m1), total_price(m1) double precision t(n1) double precision cf(m1,n1) double precision w1(m1,n1), y1(m1), A1(n1,n1), beta1(n1), b1(n1), d(n1),r(n1) double precision w2(m2,n2), y2(m2), A2(n2,n2), beta2(n2), b2(n2), a(n2) integer today, last_cpn, next_cpn integer number_of_days_since_last_cpn integer number_of_days_in_current_cpn_period ! input face_value = 1000d0 cpn_rate(1:m1)= & [6 + 5d0/8d0, 9 + 1d0/8d0, 7 + 7d0/8d0, 8 + 1d0/4d0, 8 + 1d0/4d0, & 8 + 3d0/8d0, 8d0, 8 + 3d0/4d0, 6 + 7d0/8d0, 8 + 7d0/8d0, & 6 + 7d0/8d0, 8 + 5d0/8d0, 7 + 3d0/4d0, 11 + 1d0/4d0, 8 + 1d0/2d0, & 10 + 1d0/2d0, 7 + 7d0/8d0, 8 + 7d0/8d0 ] ask_price(1:m1)= & [100+ 0d0/32d0, 100+22d0/32d0, 100+24d0/32d0, 101+1d0/32d0, 101+ 7d0/32d0, & 101+12d0/32d0, 100+26d0/32d0, 102+ 1d0/32d0, 98+5d0/32d0, 102+ 9d0/32d0, & 97+13d0/32d0, 101+23d0/32d0, 99+ 5d0/32d0, 109+4d0/32d0, 101+13d0/32d0, & 107+27d0/32d0, 99+13d0/32d0, 103+ 0d0/32d0 ] today = mjd(2011, 11, 5) last_cpn = mjd(2011, 8, 15) next_cpn = mjd(2012, 2, 15) ! calculation cpn(1:m1) = face_value * cpn_rate(1:m1) / 100d0 / 2. t(1:n1) = [(mjd(2012,2,15) - today) / 365d0, & ! unit year (mjd(2012,8,15) - today) / 365d0, & (mjd(2013,2,15) - today) / 365d0, & (mjd(2013,8,15) - today) / 365d0, & (mjd(2014,2,15) - today) / 365d0, & (mjd(2014,8,15) - today) / 365d0, & (mjd(2015,2,15) - today) / 365d0, & (mjd(2015,8,15) - today) / 365d0, & (mjd(2016,2,15) - today) / 365d0 ] number_of_days_since_last_cpn =today - last_cpn ! days number_of_days_in_current_cpn_period =next_cpn - last_cpn ! days accrued_interest(1:m1) = 1d0*(number_of_days_since_last_cpn) / & number_of_days_in_current_cpn_period * cpn(1:m1) total_price(1:m1) = face_value * ask_price(1:m1) / 100d0 + accrued_interest(1:m1) cf(1,1:n1) = [cpn( 1)+face_value,0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0] cf(2,1:n1) = [cpn( 2)+face_value,0d0,0d0,0d0,0d0,0d0,0d0,0d0,0d0] cf(3,1:n1) = [ cpn( 3),cpn( 3)+face_value,0d0,0d0,0d0,0d0,0d0,0d0,0d0 ] cf(4,1:n1) = [ cpn( 4),cpn( 4)+face_value,0d0,0d0,0d0,0d0,0d0,0d0,0d0 ] cf(5,1:n1) = [ cpn( 5),cpn( 5),cpn( 5)+face_value,0d0,0d0,0d0,0d0,0d0,0d0 ] cf(6,1:n1) = [ cpn( 6),cpn( 6),cpn( 6)+face_value,0d0,0d0,0d0,0d0,0d0,0d0 ] cf(7,1:n1) = [ cpn( 7),cpn( 7),cpn( 7),cpn( 7)+face_value, & 0d0,0d0,0d0,0d0,0d0 ] cf(8,1:n1) = [ cpn( 8),cpn( 8),cpn( 8),cpn( 8)+face_value, & 0d0,0d0,0d0,0d0,0d0 ] cf(9,1:n1) = [ cpn( 9),cpn( 9),cpn( 9),cpn( 9),cpn( 9)+ & face_value,0d0,0d0,0d0,0d0] cf(10,1:n1) = [ cpn(10),cpn(10),cpn(10),cpn(10),cpn(10)+ & face_value,0d0,0d0,0d0,0d0 ] cf(11,1:n1)= [ cpn(11),cpn(11),cpn(11),cpn(11),cpn(11), & cpn(11)+face_value,0d0,0d0,0d0 ] cf(12,1:n1)= [ cpn(12),cpn(12),cpn(12),cpn(12),cpn(12), & cpn(12)+face_value,0d0,0d0,0d0 ] cf(13,1:n1)= [ cpn(13),cpn(13),cpn(13),cpn(13),cpn(13), & cpn(13),cpn(13)+face_value,0d0,0d0 ] cf(14,1:n1)= [ cpn(14),cpn(14),cpn(14),cpn(14),cpn(14), & cpn(14),cpn(14)+face_value,0d0,0d0 ] cf(15,1:n1)= [ cpn(15),cpn(15),cpn(15),cpn(15),cpn(15), & cpn(15),cpn(15),cpn(15)+face_value,0d0 ] cf(16,1:n1)= [ cpn(16),cpn(16),cpn(16),cpn(16),cpn(16), & cpn(16),cpn(16),cpn(16)+face_value,0d0 ] cf(17,1:n1)= [ cpn(17),cpn(17),cpn(17),cpn(17),cpn(17), & cpn(17),cpn(17),cpn(17),cpn(17)+face_value ] cf(18,1:n1)= [ cpn(18),cpn(18),cpn(18),cpn(18),cpn(18), & cpn(18),cpn(18),cpn(18),cpn(18)+face_value ] w1(1:m1,1:n1)=cf(1:m1,1:n1) ! least suquare method ! (w1){beta1} = {y1} ! (w1)^T (w1){beta1} = (w1)^T {y1} ! {beta1} = ((w1)^T (w1))^-1 (w1)^T){y1} y1(1:m1) = total_price(1:m1) A1(1:n1,1:n1) = matmul(transpose(w1(1:m1,1:n1)), w1(1:m1,1:n1)) b1(1:n1) = matmul(transpose(w1(1:m1,1:n1)), y1(1:m1)) call GAUSS(A1, b1, beta1, 9, 9) d(1:n1) = beta1(1:n1) r(1:n1) = -log(d(1:n1)) / t(1:n1) ! m1xn1 -> m2xn1 do i = 1, m2 do j = 1, n2 w2(i, j) = t(i) ** (j - 1) end do end do ! least suquare method y2(1:m2) = r(1:m2) A2(1:n2,1:n2) = matmul(transpose(w2(1:m2,1:n2)), w2(1:m2,1:n2)) b2(1:n2) = matmul(transpose(w2(1:m2,1:n2)), y2(1:m2)) call GAUSS(A2, b2, beta2, 5, 5) ! output print *, "4.4.f90" print *,"face value= ", int(face_value) print *,"last_cpn= ", last_cpn print *,"today= ", today print *,"next cpn= ", next_cpn print *,"number_of_days_since_last_cpn= ", & number_of_days_since_last_cpn print *,"number_of_days_in_current_cpn_period= ", & number_of_days_in_current_cpn_period print *,"" ! print '(a5,a10,a10,a10,a15)', & "i", "a_price", "cpn", "a.int", "total_price" print *,('-', i=1,50) print '(i5,f10.2,f10.2,f10.2,f15.2)', & (i,ask_price(i),cpn(i),accrued_interest(i),total_price(i),i=1,m1) print *,('-', i=1,50) print *,"a_price = ask price" print *,"a.int = accrued interest" print *,"" ! print t print '(a10a10)', "i", "t(year)" print *,('-',i=1, 20) print '(i10,f10.3)', (i,t(i),i=1,n1) print *,('-',i=1, 20) print *,'' ! print y1 print '(a10,a10)', "i", "y1" print *,('-',i=1, 20) print '(i10,i10)', (i,int(y1(i)),i=1,m1) print *,('-',i=1, 20) print *,'' ! print w1, m1xn1-matrix print *,"w1, a matrix of cf(in integer) in row" print '(10a6)',"t0","t1", "t2","t3", "t4","t5", "t6", "t7", "t8" print *,('-',i=1,54) print '(9i6)',((int(w1(i,j)),j=1,n1),i=1,m1) print *,('-',i=1,54) print *,'' ! print beta1 print *,"estimated beta1, d(t(i))" print '(3a10)', "i", "beta1" print *,('-',i=1,20) print '(i10,f10.3)', (i,beta1(i),i=1,n1) print *,('-',i=1,20) print *,'' ! print r print *,"estimated r(i)" print '(2a10)', "i", "r(i)" print *,('-',i=1,20) print '(i10,f10.4)', (i,r(i),i=1,n1) print *,('-',i=1,20) print *,'' ! print w2, m2xn2-matrix print *,"w2, a coefficent matrix" print '(6a8)', "i","1", "t", "t^2", "t^3", "t^4" print *,("-",i=1,48) print '(i8,5f8.3)', (i,(w2(i,j),j=1,n2),i=1,m2) print *,("-",i=1,48) print *,'' ! print beta2 print *,"beta2, edtimated a(i) for r(t)" print '(a10,a15)', "i", "a(i)" print *,('-',i=1,25) print '(i10,f15.5)', (i,beta2(i),i=1,n2) print *,('-',i=1,25) print *,'' ! print r equation estimated a = beta2 print *,"estimated r(t)=" print '(f8.5, "+",f8.5," t +",f8.5," t^2 + ",f8.5," t^3 + ",f10.5," t^4")', & a(1), a(2), a(3), a(4), a(5) ! contains integer function mjd(y, m, d) ! 参照文献 ! wikipedia, ユリウス通日, グレゴリオ暦から修正ユリウス日への換算 ! https://ja.wikipedia.org/wiki/ユリウス通日 integer y, m, d , yy, mm if( m == 1) then yy = y - 1 mm = 13 else yy = y mm = m end if if(m == 2) then yy = y - 1 mm = 14 else yy = y mm = m end if mjd = floor(365.25 * yy) + floor(yy / 400.) & - floor(yy / 100.) + floor(30.59 * (mm - 2)) & + d - 678912 end end ! SUBROUTINE GAUSS(A,B,X,N,NN) ! 出典 ! 戸川隼人『数値計算』岩波書店, 1991, p.80, プログラム4.1 implicit none real(8) A(NN,NN),B(NN),X(NN),Q,S INTEGER I,J,K,N,NN DO 10 K=1,N-1 DO 20 I=K+1,N Q=A(I,K)/A(K,K) DO 30 J=K,N A(I,J)=A(I,J)-Q*A(K,J) 30 CONTINUE B(I)=B(I)-Q*B(K) 20 CONTINUE 10 CONTINUE X(N)=B(N)/A(N,N) DO 40 K=N-1,1,-1 S=B(K) DO 50 J=K+1,N S=S-A(K,J)*X(J) 50 CONTINUE X(K)=S/A(K,K) 40 CONTINUE RETURN END ! end