whittle.s

计算hurst的程序

###############################################
#
#Splus function and programs for the
#calculation of Whittle's estimator and
#the goodness of fit statistic as defi-
#ned in Beran (1992). The models are
#fractional gaussian noise or fractional
#Arima.
#The data series may be divided into
#subseries for which the parameters are
#fitted separately.
#
###############################################
#
#Functions
#
###############################################
#CetaFGN
#
#Covariance matrix of hat{eta} for fgn
###############################################

CetaFGN_function(eta)
        {
        M_length(eta)

#size of steps in Riemann sum: 2*pi/m

        m_10000
        mhalfm_trunc((m-1)/2)

#size of delta for numerical calculation of derivative

        delta_0.000000001

#partial derivatives of log f(at each fourier frequency)

        lf_matrix(1,ncol=M,nrow=mhalfm)
        f0_fspecFGN(eta,m)$fspec
        for(j in (1:M))
        {
         etaj_eta
         etaj[j]_etaj[j]+delta
         fj_fspecFGN(etaj,m)$fspec
         lf[,j]_log(fj/f0)/delta
        }

#Calculate D
        
        Djl_matrix(1,ncol=M,nrow=M)
        for(j in (1:M))
         {for(l in (1:M))
          {Djl[j,l]_2*2*pi/m*sum(lf[,j]*lf[,l])
          }
         }	

#Result

        drop(matrix(4*pi*solve(Djl),ncol=M,nrow=M,byrow=T))
        }

###############################################
#CetaARIMA
#
#Covariance matrix of hat{eta}
#for fractional ARIMA
###############################################

CetaARIMA_function(eta,p,q)
        {
        M_length(eta)

#size of steps in Riemann sum: 2*pi/m

        m_10000
        mhalfm_trunc((m-1)/2)

#size of delta for numerical calculation of derivative

        delta_0.000000001

#partial derivatives of log f (at each fourier frequency

        lf_matrix(1,ncol=M,nrow=mhalfm)
        f0_fspecARIMA(eta,p,q,m)$fspec
        for(j in (1:M))
         {
         etaj_eta
         etaj[j]_etaj[j]+delta
         fj_fspecARIMA(etaj,p,q,m)$fspec
         lf[,j]_log(fj/f0)/delta
         }

#Calculate D

        Djl_matrix(1,ncol=M,nrow=M)
        for(j in (1:M))
         {for(l in (1:M))
          {Djl[j,l]_2*2*pi/m*sum(lf[,j]*lf[,l])
          }
         }

#Result

        drop(matrix(4*pi*solve(Djl),ncol=M,nrow=M,byrow=T))
        }

###############################################
#Qeta
#
#Function for the calculation of A, B and
#Tn=A/B**2
#where A = 2pi/n sum 2*[I(lambda_j)/f(lambda_j)],
#B = 2pi/n sum 2*[I(lambda_j)/f(lambda_j)]**2 ans
#the sum is taken over all fourier frequencies
#lambda_j=2pi*j/n (j=1,...,(n-1)/2.
#f is the spectral density of fractional gaussian
#noise or fractional ARIMA(p,d,q)
#with self-similarity parameter H=h.
#cov(X(t),X(t+k))=integral(exp(iuk)f(u)du)
#
#NOTE: yper[1] must be the periodogram I(lambda_1) at
#the frequency @pi/n (i.e. not the frequency zero!).
#
#INPUT: h
#
#(n,nhalfm=trunc[(n-1)/2] and the
#nhalfm-dimensional vector yper must
#be defined.)
#
#
#OUTPUT: list(n=n,h=h,A=A,B=B,Tn=Tn,z=z,pval=pval,
#theta1=theta1,fspec=fspec)
#
#Tn is the goodness of fit test statistic
#Tn=A/B**2 defined in Beran (1992),
#z is the standardized test statistic,
#pval the corrsponding p-value P(w>z).
#theta1 is the scale parameter such that
#f=theta1*fspec and integral(log[fspec])=0
#
###############################################

Qeta_function(eta)
        {
#       cat("in function Qeta", fill=T)
        h_eta[1]

#spectrum at fourier frequencies

        if(imodel==1)
         {fspec_fspecFGN(eta,n)
         theta1_fspec$theta1
         fspec_fspec$fspec}
        else
         {fspec_fspecARIMA(eta,p,q,n)
         theta1_fspec$theta1
         fspec_fspec$fspec}

#Tn=A/B**2

        yf_yper/fspec
        yfyf_yf**2
        A_2*(2*pi/n)*sum(yfyf)
        B_2*(2*pi/n)*sum(yf)
        Tn_A/(B**2)
        z_sqrt(n)*(pi*Tn-1)/sqrt(2)
        pval_1-pnorm(z)
        theta1_B/(2*pi)
        fspec_fspec
        Qresult_list(n=n,h=h,
                eta=eta,A=A,B=B,Tn=Tn,z=z,pval=pval,
                theta1=theta1,fspec=fspec)
        drop(Qresult)
        }

###############################################
#calculation of the spectral density of
#normalized fractional gaussian noise
#with self-similarity parameter H=h
#at the Fourier frequencies 2*pi*j/m (j=1,...,(m-1)).
#
#Remarks:
#
#1. cov(X(t),X(t+k))=integral[exp(iuk)f(u)du]
#2. f=theta1*fspec and integral[log(fspec)]=0
#
#INPUT: m=sample size
#h=self-similarity parameter
#
#OUTPUT:list(fspec=fspec,theta1=theta1)
#
###############################################








fspecFGN_function(eta,m)
        {
#parameters for the calculation of f       
#        h_eta[1]
#        nsum_200
#        hh_-2*h-1
#        const_1/pi*sin(pi*h)*gamma(-hh)
#        j_2*pi*c(0:nsum)

#x=2*pi*(j-1)/m (j=1,2,...,(n-1)/2)
#Fourier frequency

#        mhalfm_trunc((m-1)/2)
#        x_1:mhalfm
#        x_2*pi/m*x

#calculation of f at fourier frequencies

#        fspec_matrix(0,mhalfm)
#        for(i in seq(1:mhalfm))
#         {
#         lambda_x[i]
#         fi_matrix(lambda,(nsum+1))
#         fi_abs(j+fi)**hh+abs(j-fi)**hh
#         fi[1]_fi[1]/2
#         fi_(1-cos(lambda))*const*fi
#         fspec[i]_sum(fi)
#         }
	fspec <- ffourier.FGN.est(eta, m)
#Above change made to speed up routine.  Due to Vern Paxson.   3/23/95 VT
###############################################

        logfspec_log(fspec)
        fint_2/(m)*sum(logfspec)
        theta1_exp(fint)
        fspec_fspec/theta1
        drop(list(fspec=fspec,theta1=theta1))
        }

ffourier.FGN.est <-function(H, n)
FGN.spectrum((2 * pi * (1:((n - 1)/2)))/n, H)/pi/2

FGN.spectrum <-function(lambda, H)
{
        2 * sin(pi * H) * gamma(2 * H + 1) * (1 - cos(lambda)) * (lambda^(-2 * 
                H - 1) + FGN.B.est.adjust(lambda, H))
}

FGN.B.est.adjust <- function(lambda, H)
{
        B <- FGN.B.est(lambda, H)
        (1.0002 - 0.000134 * lambda) * (B - 2^(-7.65 * H - 7.4))
}

FGN.B.est <- function(lambda, H)
{
        d <-  - (2 * H + 1)
        dprime <- -2 * H
        a <- function(lambda, k)
        2 * k * pi + lambda
        b <- function(lambda, k)
        2 * k * pi - lambda
        a1 <- a(lambda, 1)
        b1 <- b(lambda, 1)
        a2 <- a(lambda, 2)
        b2 <- b(lambda, 2)
        a3 <- a(lambda, 3)
        b3 <- b(lambda, 3)
        a4 <- a(lambda, 4)
        b4 <- b(lambda, 4)
        a1^d + b1^d + a2^d + b2^d + a3^d + b3^d + (a3^dprime + b3^dprime + a4^
                dprime + b4^dprime)/(8 * pi * H)
}


###############################################
#calculation of the spectral density of
#fractional ARMA with standard normal innovations
#and self similarity parameter H=h
#at the fourier frequencies 2*pi*j/n (j=1,...,(n-1)).
#cov(X(t),X(t+k))=(sigma/(2*pi))*integral(exp(iuk)g(u)du).
#
#INPUT: m=sample size
#h=theta[1]=self-similarity parameter
#phi=theta[2:(p+1)]=AR(p)-parameters
#psi=theta[(p+2):(p+q+1)]=MA(q)-parameters
#
#OUTPUT:list(fspec=fspec,theta1=theta1)
#
###############################################

fspecARIMA_function(eta,p,q,m)
        {
###############################################
#       cat("in fspecARIMA",fill=T)
        h_eta[1]
        phi_c()
        psi_c()

#x=2*pi*(j-1)/m (j=1,2,...,(n-1)/2)
#=Fourier frequencies

        mhalfm_trunc((m-1)/2)
        x_1:mhalfm
        x_2*pi/m*x

#calculation of f at fourier frequencies

        far_(1:mhalfm)/(1:mhalfm)
        fma_(1:mhalfm)/(1:mhalfm)
        if(p>0)
         {phi_cbind(eta[2:(p+1)])
         cosar_cos(cbind(x)%*%rbind(1:p))
         sinar_sin(cbind(x)%*%rbind(1:p))
         Rar_cosar%*%phi
         Iar_sinar%*%phi
         far_(1-Rar)**2+Iar**2
         }
#       cat("far calculated",fill=T)

        if(q>0)
         {psi_cbind(eta[(p+2):(p+q+1)])
         cosar_cos(cbind(x)%*%rbind(1:q))
         sinar_sin(cbind(x)%*%rbind(1:q))
         Rar_cosar%*%psi
         Iar_sinar%*%psi
         fma_(1+Rar)**2+Iar**2
         }
#       cat("fma calculated", fill=T)

        fspec_fma/far*sqrt((1-cos(x))**2+sin(x)**2)**(1-2*h)
        theta1_1/(2*pi)
#       cat("end of fspecARIMA",fill=T)
        drop(list(fspec=fspec,theta1=theta1))
        }

###############################################
#definition of the periodogram
###############################################

per_function(z)
        {n_length(z)
        (Mod(fft(z))**2/(2*pi*n)) [1:(n%/%2+1)]
        }

###############################################
#definition of function to be minimized
###############################################

Qmin_function(etatry)
        {
        result_Qeta(etatry)$B
	if(out==T)
	{
        cat("etatry=",etatry,"B=",result,sep=" ",fill=T)}
	assign("bBb",result,frame=0)
        drop(result)
        }

###############################################
#           M A I N    P R O G R A M
###############################################

whittle_function(xinput, nsub = 1, model = "farima", pp = 1, qq = 1, h =
0.5, ar = c(0.5), ma = c(0.5), out = T, spec = F)
{
	assign("out", out, frame = 1)
	nmax <- length(xinput)
	startend <- c(1, nmax)
	istart <- startend[1]
	iend <- startend[2]
	nloop <- nsub
	assign("n", trunc((iend - istart + 1)/nloop), frame = 1)
	nhalfm <- trunc((n - 1)/2)
	if(model == "farima")
		assign("imodel", 2, frame = 1)
	if(model == "fgn")
		assign("imodel", 1, frame = 1)
	assign("p", 0, frame = 1)
	assign("q", 0, frame = 1)
	if(imodel == 2) {
		assign("p", pp, frame = 1)
		assign("q", qq, frame = 1)
	}
	else {
		assign("p", 0, frame = 1)
		assign("q", 0, frame = 1)
	}
	eta <- c(h)
	if(p > 0) {
		eta[2:(p + 1)] <- ar
	}
	if(q > 0) {
		eta[(p + 2):(p + q + 1)] <- ma
	}
	M <- length(eta)	#loop
	thetavector <- c()
	i0 <- istart
	flax <- vector("list", nloop)	

#necessary to make nsub/nloop work.  VT.

	for(iloop in (1:nloop)) {
		h <- max(0.2, min(h, 0.9))
		eta[1] <- h
		i1 <- i0 + n - 1
		y <- xinput[i0:i1]	

#standardize data

		vary <- var(y)
		y <- (y - mean(y))/sqrt(var(y))	

#periodogram of the data

		if(spec == T) {
			assign("yper", per(y)[2:(nhalfm + 1)], w = 1)
		}
		else {
			assign("yper", per(y)[2:(nhalfm + 1)], frame = 1)
		}
		s <- 2 * (1 - h)
		etatry <- eta	

#	Modified to make nlmin not give incorrect result.  VT

		if(imodel == 1) {
			result <- nlminb(start = etatry, objective = Qmin, 
				lower = 0, upper = 0.999)
		}
		else {
			result <- nlminb(start = etatry, objective = Qmin)
		}

		eta <- result$parameters
		sturno <- result$message
		theta1 <- Qeta(eta)$theta1
		theta <- c(theta1, eta)
		thetavector <- c(thetavector, theta)	

#calculate goodness of fit statistic

		Qresult <- Qeta(eta)	#output
		M <- length(eta)
		if(imodel == 1) {
			SD <- CetaFGN(eta)
			SD <- matrix(SD, ncol = M, nrow = M, byrow = T)/n
		}
		else {

#  Changed to eliminate crashing in solve.qr in CetaARIMA.  VT

			cat("M = ", M, "\n")
			if(M > 2) {
				for(i in 3:M) {
				  for(j in 2:(i - 1)) {
				    temp <- eta[i] + eta[j]
				    if(abs(temp) < 0.0001) {
				      cat(
	        "Problem with estimating confidence intervals, parameter ",
				        i, "and  parameter ", j, 
				        "are the same, eliminating.\n")
				      eta <- eta[ - i]
				      eta <- eta[ - j]
				      M <- M - 2
				      p <- p - 1
				      q <- q - 1
				    }
				  }
				}
			}
			SD <- CetaARIMA(eta, p, q)
			SD <- matrix(SD, ncol = M, nrow = M, byrow = T)/n
		}
		Hlow <- eta[1] - 1.96 * sqrt(SD[1, 1])
		Hup <- eta[1] + 1.96 * sqrt(SD[1, 1])
		if(out == T) {
			cat("theta=", theta, fill = T)
			cat("H=", eta[1], fill = T)
			cat("95%-C.I. for H: [", Hlow, ",", Hup, "]", fill = T)
		}
#
#Changing of the signs of the moving average parameters
#in order to respect the sign of the Splus convention
#
		if(q > 0) {
			eta[(p + 2):(p + q + 1)] <-  - eta[(p + 2):(p + q + 1)]
		}
		etalow <- c()
		etaup <- c()
		for(i in (1:length(eta))) {
			etalow <- c(etalow, eta[i] - 1.96 * sqrt(SD[i, i]))
			etaup <- c(etaup, eta[i] + 1.96 * sqrt(SD[i, i]))
		}
		if(out == T) {
			cat("95%-C.I.:", fill = T)
			print(cbind(etalow, etaup), fill = T)
		}
		if(spec == T) {
			cat("periodogram is in yper", fill = T)
			assign("fest", Qresult$theta1 * Qresult$fspec, w = 1)
			cat("spectral density is in fest", fill = T)
		}
		flax[[iloop]] <- list()		# To make nsub work.  VT.
		flax[[iloop]]$parameter <- eta
		flax[[iloop]]$sigma2 <- bBb * var(xinput)
		flax[[iloop]]$conv.type <- sturno
		remove("bBb", frame = 0)		#next subseries
		i0 <- i0 + n	
#
#Changing of the signs of the moving average parameters
#
		if(q > 0) {
			eta[(p + 2):(p + q + 1)] <-  - eta[(p + 2):(p + q + 1)]
		}
	}
	return(flax)
#
#	Clean up
#
	remove("n", frame = 1)
	remove("p", frame = 1)
	remove("q", frame = 1)
	remove("imodel", frame = 1)
	remove("out", frame = 1)
	if(spec == F) {
		remove("yper", frame = 1)
	}
}