# plotmoments.S
#Written by Vadim Teverovsky
	
# These functions use the absolute moments of a series to estimate H.
# plotmoment subtracts the mean and calculates the moments.
# plotmoment2 computes the moments without first subtracting the mean.  If a
# series has non-zero mean, the second method may give incorrect estimate, and
# non-linear plot. 

plotmoment_function(data, nvar = 50, minnpts = 5, power1 = 0.7, power2 =
2.5, slopes = 1,  plotflag = T, moment = 2)
{

# This is a function which gives a plot of the moments(absolute) for
# various aggregation levels.  It calls the function Smoments.S, which,
# in turn, calls Cmoments.c.  The estimate of the slope, beta, gives
# an estimate for the parameter, H, via :  H = (beta+moment)/moment. 
#
# Variables:
#         data    --      The data being used.
#         nvar    --      The number of aggregation levels.
#         minnpts --      The minimum number of points to be used to
#                         estimate the variance at any aggregation level
#         power1,power2-- Cut-offs for points to be used in estimating
#                         beta.
#         slopes  --      Flag.  If 1, draw the estimated line whose 
#                         slope is beta.
#         plotflag  --    Flag.  If T, plot is performed.  12/27/94
#         moment   --     What absolute moment to calculate.

        cat("\n")

# n is the length of the data vector, variances is the vector of length nvar+1
# which is the output of Smoments.  It is the vector of moments at the 
# different levels of aggregation.

        n <- length(data)
        variances <- Smoments(data, nvar, minnpts, moment)

# The increment to be used in plotting and estimating the slope.  Depends on 
# the number of levels, and the minimum number of points needed at each level.

        increment <- (log10(n/minnpts))/nvar

# lv and lm are going to be the y and x values respectively of the points in
# the plot.

        lv <- log10(variances)
        lm <- rep(0, nvar + 1)

# power1 and power2 determine whether the points will be included in the 
# estimate of beta.

        for (i in 1:(nvar + 1))
          lm[i] <- log10(floor(10^((i - 1) * increment)))

# Above line changed 2/28/95, to make the plotting consistent with
# calculations. 

        lmadj <- lm[(lm >= power1) & (lm <= power2)]
        lvadj <- lv[(lm >= power1) & (lm <= power2)]
        lmadjc <- lm[(lm < power1) | (lm > power2)]
        lvadjc <- lv[(lm < power1) | (lm > power2)]

        # The line is fitted to the appropriate points.

        fit <- lsfit(lmadj, lvadj)
        intercept <- fit$coef[1]
        beta <- fit$coef[2] # The plot is done here.

# Changed the normalization from dividing by lv[1] to subtracting it.(Logs!)
# 5/23/1994.  Plotflag added 12/27/94

        if (plotflag == T)
          {
            plot(c(0, lm), c(0, lv - lv[1]), xlab = "log10(m)", ylab = 
                 "log10(moments)", type = "n")
            points(lmadj, lvadj - lv[1], pch = 16, mkh = 0.06)
            points(lmadjc, lvadjc - lv[1], pch = 3, mkh = 0.06)
        #          abline(0,  - moment/2, lty = 2)
            if (slopes == 1)
              {
                lines(lmadj, intercept - lv[1] + (beta) * lmadj)
              }
        }

# Returning results 12/27/94.

        cat(" beta = ", beta, "\n\n")
        cat("H= ", (beta + moment)/moment, "\n\n")
        result <- (beta + moment)/moment
        names(result) <- NULL
        return(result, variances)
}

Smoments
function(data, nvar, minpts, moment)
{
        n <- length(data)
        if(!is.loaded(C.symbol("Cmoments")))
                dyn.load("Cmoments.o")
        .C("Cmoments",
                data = as.double(data),
                n = n,
                nvar = as.integer(nvar),
                minpts = as.integer(minpts),
                moment = as.double(moment),
                v = double(nvar + 1))$v
}

plotmoment2 _ 
function(data, nvar = 50, minnpts = 5, power1 = 0.7, power2 = 2.5, slopes = 1, 
        plotflag = T, moment = 2)
{

#       This is a function which gives a plot of the variances for various
#       aggregation levels.  It calls the function Smoments2.S, which, in
#       turn, calls Cmoments2.c.  The estimate of the slope, beta, gives
#       an estimate for the parameter, H, via :  H = (beta+2)/2.
#       Variables:
#
#               data    --      The data being used.
#               nvar    --      The number of aggregation levels.
#               minnpts --      The minimum number of points to be used to
#                               estimate the variance at any aggregation level
#               power1,power2-- Cut-offs for points to be used in estimating
#                               beta.
#               slopes  --      Flag.  If 1, draw the estimated line whose 
#                               slope is beta.
#               plotflag  --    Flag.  If T, plot is performed.  12/27/94
#

        cat("\n")

# n is the length of the data vector, variances is the vector of length nvar+1
# which is the output of Svariances.  It is the vector of variances at the 
# different levels of aggregation.

        n <- length(data)
        variances <- Smoments2(data, nvar, minnpts, moment)
# The increment to be used in plotting and estimating the slope.  Depends on 
# the number of levels, and the minimum number of points needed at each level.

        increment <- (log10(n/minnpts))/nvar

# lv and lm are going to be the y and x values respectively of the points in
# the plot.

        lv <- log10(variances)
        lm <- rep(0, nvar + 1)

# power1 and power2 determine whether the points will be included in the 
# estimate of beta.

        for (i in 1:nvar + 1)
                lm[i] <- log10(floor(10^((i - 1) * increment)))
# Above line changed 2/28/95, to make the plotting consistent with
# calculations.

        lmadj <- lm[(lm >= power1) & (lm <= power2)]
        lvadj <- lv[(lm >= power1) & (lm <= power2)]
        lmadjc <- lm[(lm < power1) | (lm > power2)]
        lvadjc <- lv[(lm < power1) | (lm > power2)]

        # The line is fitted to the appropriate points.

        fit <- lsfit(lmadj, lvadj)
        intercept <- fit$coef[1]
        beta <- fit$coef[2] # The plot is done here.

# Changed the normalization from dividing by lv[1] to subtracting it.(Logs!)
# 5/23/1994.  Plotflag added 12/27/94

        if (plotflag == T)
          {
            plot(c(0, lm), c(0, lv - lv[1]), xlab = "log10(m)", ylab = 
                 "log10(moments)", type = "n")
            points(lmadj, lvadj - lv[1], pch = 16, mkh = 0.06)
            points(lmadjc, lvadjc - lv[1], pch = 3, mkh = 0.06)
        #                  abline(0,  - moment/2, lty = 2)
            if (slopes == 1)
	      {
                lines(lmadj, intercept - lv[1] + (beta) * lmadj)
              }
          }

# Returning results 12/27/94.

        cat(" beta = ", beta, "\n\n")
        cat("H= ", (beta + moment)/moment, "\n\n")
        result <- (beta + moment)/moment
        names(result) <- NULL
        return(result)
}

Smoments2
function(data, nvar, minpts, moment)
{
        n <- length(data)
        if(!is.loaded(C.symbol("Cmoments2")))
                dyn.load("Cmoments2.o")
        .C("Cmoments2",
                data = as.double(data),
                n = n,
                nvar = as.integer(nvar),
                minpts = as.integer(minpts),
                moment = as.double(moment),
                v = double(nvar + 1))$v
}

# eof plotmoments.S
