S-PLUS

The statistical methods described here use S-Plus, a statistical package available at http://www.mathsoft.com/splus.html.

Additional software and extensions for S-Plus are available at http://lib.stat.cmu.edu/S/.

The command S-Plus will start up the package on most platforms. This will automatically create a .Data directory where the S-Plus data will be stored. The examples shown in these pages assume that whatever files are being scanned in are in the same directory where S-Plus was started, as are any object files which are necessary. If this is not the case, the path should be included in the scan calls. The C files which are used in some of the methods above need to be compiled to an object file, which is what the S-Plus routines actually use.

Splus is a statistical package which allows the user a wide flexibility in
usage.  

Simple Splus objects include not only numbers, vectors and matrices, but
also inherent and user-defined functions.

Following is a sample session, with comments denoted by "#". ">" is the prompt.
Splus	# Starts Splus from shell

X11()   

#Starts graphics window.

#All of the usual arithmetic and logical operators are present.  
5+3
> 8

#The "help()" function allows for on-line help about any Splus object, be it
#function, etc.

help(X11) #Provides help on the command X11().

#To attach a directory different form .Data, use the attach command:
attach("/usr/joe/.Data")

#For more information
help(attach)


#One can create and dimension a vector of length n, filled with zeros, by
#the command "vect <- rep(0,n)".   "<-" or "_" is the assignment operator.  

hestimates <- rep(0,50)

#This is a nice way to create a vector which can be used later, e.g. in a
#FOR loop.  The vector elements are referred to as "hestimates[i]".


#One can simulate stable series with the command "rstab".

sample.series <- rstab(1000, index=1.5)  # Length 1000, alpha=1.5.

tsplot(sample.series)	# Produces time series plot.

#One can also simulate non-symmetric stable series -- the parametrization 
#used in this case is special -- see Samorodnitsky and Taqqu, Chapman 
#and Hall, 1994, p. 43.

sample.series[1]
> 0.8194

lsfit(sample.series)$coefficients  
> 1   .5  
# $coefficients produces a subset of the output, namely, the coefficients
# of a least square fit(Intercept, slope). 
# (Residuals, etc. can also be obtained).

stem(sample.series)  # Produces stem plot (median, quartiles, etc.)

# One can simulate long-range dependent series (FARIMA) with the command
# "arima.fracdiff.sim" .  For the right syntax, use help:

help(arima.fracdiff.sim)

#or look at the function:

arima.fracdiff.sim

sample.series2 <- arima.fracdiff.sim(model=list(d=0.3),n=1000)
sample.series2[1]
> 1.293

#Most of the statistical methods take a vector as the "data" variable, and
#run the vector through some manipulations, and usually come up with
#numerical results and a plot.  For example, 

#(1) Create and dimension an array:
farima.data_array(0,dim=c(50,1000))  
# Creates an array of zeros 50 by 1000.  Members of array can be referenced as 
# farima.data[i,j].  A whole column or row of array can be referenced as 
# farima.data[i,] or  farima.data[,j].

#(2) Fill in array 
# Now we put simulated long-range dependent series into array, making 50
# series, each length 1000.  We use a FOR loop.

for(i in 1:50){
	 farima.data[i,]<-arima.fracdiff.sim(model=list(d=0.1),n=1000)}

#(3) Use an estimation procedure:
# Suppose we want to estimate the scaling parameter H using the Aggregated
# Variance method.  We will use the function "plotvar" contained in the
# file "file.S".  The source command reads in the contents of this file.

source("file.S")

#Now we estimate H using our Aggregated Variance method.
for (i in 1:50}{
	hestimates[i]_plotvar(farima.data[i,],plot=False)}

#(4) Examine results:
#Since we are analyzing 50 series, we do not usually want to see plots.

#Now one can find statistics for the estimator of H.
mean(hestimates)
sqrt(var(hestimates))
hist(hestimates) #returns histogram of this vector.

#(5) Writing results to outside file, sink all output:
sink("file_name")
mean(hestimates)
sqrt(var(hestimates))
sink() # closes file.


# To see what objects one has defined, use the command:

objects()  

# This lists everything which has been defined and not deleted.
# To delete things use rm().

rm(sample.series)


# To see what an object is, type in the name:

plotvar	#returns the definition of the function plotvar
hestimates  #returns the values of the vector hestimates


#To split graphics window into 4 sub-plots, 2 by 2.  

par(mfrow=c(2,2))

#For info on axes, title, captions, etc.

help(plot)

#To save a graphic (plot): One has to regenerate the output and direct it
#to a file instead of the X11 window.  Last device created and not yet
#closed is the default.   Open file:

postscript("name_of_file.ps")
plotvar(sample.series, plot=T)  #generate plot  (portrait is default).
dev.off()                       #close file 

#Alternatively, the print button on the X11() window saves a postscript
#file as a by-product of printing.  (Saved in directory where S-Plus was 
#started.  The default is landscape, but portrait printing can be done by 
#pressing the "print orientation" button in the X11() window.

# Control-C interrupts.  Control-Z suspends.

#Quit

q()