dfa_ny.f90

给出去趋势分析DFA方法的精确求解

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!          This program is for Detrended Fluctuation Analysis           !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!




!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
! parameter nfit denotes the order of polynominals used in the process of fitting
! minbox: minimum boxsize; maxbox: maximum boxsize
! npts: the length of the original time series 
!boxratio: the ratio between the sequential boxes
!Ys():the difference between the original time series and the fits in the each box
!Fs(:): the variance--Fs()^2=<Ys()^2>
!F(): root-mean-square (rms) fluctuation for the detrended signals, named, DFA fluctuation funtion
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
program main
use IMSL
implicit none
integer,parameter :: nfit=2
real(kind=8) :: boxratio,buffer
integer :: minbox, maxbox,i,j,k,m,v,n,rslen,boxsize,npts,ERROR,count,status=0,error1 
real(kind=8) :: coe(nfit+1),b0,tempy,temp,xaver(nfit),poly(nfit+1), stat(12)
character(len=5) :: name
character(len=30) :: inname
real(kind=8),allocatable :: xx(:),yy(:),Fs(:),F(:),seq(:), Y(:),pv(:),Ys(:)


! file station.txt contains the information of 194 stations
open(888,file="f:\study\allstations191.txt")

do while(.true.)
read(888,"(A5)",iostat=error1) name

!  if(error1/=0)exiteal(kind=8),
  inname="f:\study\rizhaoaver\s"//name//"_dev"
  print*,inname
! file 99 save the output data log10(F(n)),log10(real(boxsize))
open(unit=99,file=("f:\study\out\out"//name//"_ny2.dat"))
open(777,file="f:\study\profilerh\profilerh"//name//".dat")
! file 10 input the original data
open(unit=10,file=inname//".dat", iostat=ERROR)
!open(unit=10,file="f:\study\rizhaoaver\s50353_dev.dat", iostat=ERROR)
if (ERROR/=0)then	 !check whether the file exists
  write(*,*)"open file fail!"
  stop
end if
! count the the length of the time series inputted
count=0
do while(.true.)
  read (10,100,iostat=status) buffer
  100 format(F10.3) ! set the reading format
  if(status/=0)	     exit
  if((buffer/=32766).AND.(buffer/=32744)) then
  count=count+1
  end if
end do
npts=count
print*,npts
minbox=8*nfit
maxbox=npts/4

rewind(10)
allocate(seq(npts))
allocate(Y(npts))
allocate(pv(npts))
allocate(Ys(npts))
Y(:)=0.0
count=1
! put the original time series to the array: seq()
do while(.true.)
  read (10,100,iostat=status) buffer
  if(status/=0)	     exit
  if(buffer/=32766.and.buffer/=32744) then
  seq(count)=buffer
  count=count+1
  end if
end do
! subroutine deter() to compute the profile Y(i)=sum(seq(1:i))-<seq(:)>
call deter(seq,npts,Y)
do i=1,npts
  write(777,"(1x,f13.5)") Y(i)
end do
close(777)
boxratio=2.0**(1.0/5.0)
boxsize=minbox
n=1
! estimate the length of F() array
print*,name

rslen=int(log10(real(maxbox/minbox))/log10((boxratio))+10.5)


write(*,*) rslen
allocate(F(rslen))
! for each boxsize, compute the F()
do while(boxsize<=maxbox)
  v=1 
  allocate(Fs(2*int(npts/boxsize)))
  allocate(yy(boxsize))
  allocate(xx(boxsize))
! compute the F() from 1 to npts
  do i=1,npts,boxsize
 	 if(i+boxsize>npts+1) exit
	   do j=1,boxsize
	      xx(j)=real(i+j-1)
	   end do 
	 yy(1:boxsize)=Y(i:(i+boxsize-1))
    ! use the least-square method to fit the profile Y()
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  subroutine drcurv: to fit a polynomial curve using least squares.                           !
!  Arguments:                                                                                  !
!  boxsize: Number of observations. (Input)                                                    !
!  xx: Vector of length boxsize containing the x values. (Input)                               !
!  yy: Vector of length boxsize containing the y values. (Input)                               !
!  nfit: Degree of polynomial. (Input)                                                         !
!  coe: Vector of length boxsize + 1 containing the coefficients. (Output)                     !
!  poly: Vector of length boxsize + 1 containing the sequential sums of squares. (Output)      !
!  stat: Vector of length 10 containing statistics described below. (Output)                   !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
	 call DRCURV(boxsize,xx,yy,nfit,coe,poly,stat)
     Fs(v)=0.0
	 do m=i,i+boxsize-1
	   pv(m)=coe(1)
	   do j=1,nfit
          temp=1.0	  
		  do k=1,j
		     temp=temp*m
          end do
		 pv(m)= pv(m)+ coe(j+1)*temp
	   end do	
	   Ys(m)=Y(m)-pv(m)
	   Fs(v)=Fs(v)+Ys(m)*Ys(m)
	 end do	
	 Fs(v)=Fs(v)/real(boxsize)
	 v=v+1		 
   end do	
  ! compute the F() again, from npts to 1
   do i=npts,1,-boxsize
     if(i-boxsize<0) exit
	 do j=1,boxsize
	   xx(j)=real(i-j+1)
	 end do	 
	 do j=1,boxsize
	    yy(j)=Y(i-j+1)
	 end do
	 tempy=0.0
	 call DRCURV(boxsize,xx,yy,nfit,coe,poly,stat)
	 Fs(v)=0.0
	 do m=i,i-boxsize+1,-1
	   pv(m)=coe(1)
	   do j=1,nfit
	     temp=1	  
		  do k=1,j
		     temp=temp*m
		  end do   
	     pv(m)= pv(m)+ coe(j+1)*temp
	   end do	
	   Ys(m)=Y(m)-pv(m)
	   Fs(v)=Fs(v)+Ys(m)*Ys(m)
	 end do
	 Fs(v)=Fs(v)/real(boxsize)
	 v=v+1
  end do
  F(n)=0.0
  do i=1,2*int(npts/boxsize)
  F(n)=F(n)+Fs(i)
  end do
  F(n)=(F(n)/real(2*int(npts/boxsize)))
!  write(99,"(F11.6,F11.6)") 1.0/2.0*log10(F(n)),log10(real(boxsize))
 write(99,"(F11.6,F11.6)") log10(real(boxsize)),1.0/2.0*log10(F(n))
  ! write(*,*) log10(F(n)), log10(real(boxsize))
  deallocate(xx)
  deallocate(yy)
  deallocate(Fs)
  boxsize=int(boxsize*boxratio)
  n=n+1
end do
write(*,*) "rslen=", n
close(99)
close(10)  
deallocate(F)
deallocate(seq)
deallocate(Y)
deallocate(pv)
deallocate(Ys)
end do
close(888)
end program main
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! subroutine deter to  compute the profile Y(i)=sum(seq(1:i))-<seq(:)>!     !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine deter(seq,n,Y) 	
implicit none
integer :: i,j,n
real(kind=8) :: seq(n),aver,sum,Y(n)
sum=0.0
do i=1,n
  sum=sum+seq(i)
end do
aver=sum/real(n)
do i=1,n
    Y(i)=0.0
  do j=1,i
    Y(i)=Y(i)+(seq(j)-aver)
  end do
end do
return
end	subroutine deter