farima.generate <- function(n, N, d, alpha, theta, phi, sigma, beta = 0)
{
# This is a function which generates FARIMA(1,d,1) processes with stable
# innovations.  It can easily # be modified to generate any other
# FARIMA(p,d,q) processes by modifying the function armasim.  
# 
# This function calls:  
#	rstab  -- S-Plus routine which generates i.i.d. stable random
#		variables. 
#	armasim -- function which generates an ARMA series.  It calls:
#		first.x  -- function which generates the first term in series.
#	kernel.gen -- function which generates the coefficients in the
#		expansion of the fractional difference.  This is
#		independent of  everything except n and d.
#	first.x -- function to generate the right variance for the first term.
# n is the length of the kernel, or the length of summation, in X_i =
# sum_{j=0}^(n-1).  N is the legnth of desired series.   TT2 is length of
# noise vector which has to be simulated.  d is fractional difference
# parameter.  alpha is index, beta is the skewness of the stable
# distribution, theta and phi are MA and AR coefficients respectively.  sigma
# is the scale parameter of the innovations. 

# find length of needed noise vector

        TT2 <- n + N - 1

# if no MA and AR components, just do the if, if there are, the else.

        if(theta == 0 & phi == 0) {
                e <- sigma * (rstab(TT2, alpha, beta) + beta*tan(pi*alpha/2))
        }
        else {
                e <- armasim(TT2, alpha = alpha, beta = beta, ar = phi, ma = 
                        theta, sigma = sigma)
        }


        k <- kernel.gen(n, d)
# Fractional differencing.

        x <- filter(e, k, method = "convolution", sides = 1)

        y <- x[n:TT2]
        return(as.vector(y))
}

 armasim <- function(n, ar, ma, alpha, beta, sigma)
{
# create innovations

        z <- sigma * (rstab(n, index = alpha, skew = beta) +
						beta*tan(pi*alpha/2)) 
        x <- rep(0, n)

# generate first term in series

        x[1] <- (rstab(1, index = alpha, skew =beta)+beta*tan(pi*alpha/2))* first.x(a = alpha, s = sigma, th = ma, phi =ar) 
# generate the rest
        for(i in 2:n)
                x[i] <- z[i] - z[i - 1] * ma + ar * x[i - 1]
        return(x)
}

first.x <- function(alpha = 2, sigma = 1, phi = 0, theta = 0)
{
        x <- sigma * (1 + (abs(phi - theta))^alpha/(1 - (abs(phi))^alpha))^(1/
                alpha)
        return(x)
}

kernel.gen  <- function(n, d)
{
#This function generates the kernel of desired length, n, and with the
#desired long-range parameter d.  This kernel is the expansion of
#delta^(-d), which is used to generate a FARIMA(p,d,q) sequence  from a 
#ARMA(p,q,) sequence using a convolution.  See Samorodnitsky and Taqqu,
#1994, p.381. 

        b <- rep(0, (n - 1))
        k <- 1:(n - 1)
        k1 <- (k + d - 1)/k
        for(i in 1:(n - 1))
                b[i] <- prod(k1[1:i])
        b <- c(1, b)
        return(b)
}


farima.generate.pareto<-function(n, N, d, alpha, sigma)
{
# Thiis a function which generates FARIMA(0,d,0) processes. Pareto innovations.
# 
# This function calls:  
#	kernel.gen -- function which generates the coefficients in the
#		expansion of the fractional difference.  This is
#		independent of  everything except n and d.
# n is the length of the kernel, or the length of summation, in X_i =
# sum_{j=0}^(n-1).  N is the legnth of desired series.   TT2 is length of
# noise vector which has to be simulated.  d is fractional difference
# parameter.  alpha is index, sigma is the scale parameter of the innovations.
# find length of needed noise vector
	TT2 <- n + N - 1
	e <- sigma * runif(TT2)^(-1/alpha)
	k <- kernel.gen(n, d)	# Fractional differencing.
	x <- filter(e, k, method = "convolution", sides = 1)
	y <- x[n:TT2]
	return(as.vector(y))
}