This directory contains the S code (together with the C code called by
the S functions) and the original data used in the paper:

R.A. Carmona and A. Wang:
Comparison tests for the spectra of dependint multivariate time series
  
which appeared in Stochastic Modelling in Physical Oceanography, eds
R.J. Adler, P. Muller and B.L. Rozovskii, Birkhauser (1996).

It also contains the postscript files included as figures in the paper
as well as the numerical results of the tests reproduced in the tables
of the paper.


------------------------------------------------------------------------------

                  Description of the S functions

------------------------------------------------------------------------------

file name		S functions
---------		-----------
read_vel.s		make.vel.frame()
			read.vel()

scan_drifter.s		scan.drifter( filename )

plot_drifter.s		plot.drifter( drifter.list, scale=F )
			plot.inter( ... )
			drifter.world()
			motion( starting=2446200, n=100, plot.all=F )

wishart.s		wishart.ratio( n=100, r=2, K=10, method="determinant" )
			det( x )
			wishart.sample()
			myqq( sample1, sample2, 
				xlab=deparse(substitute(sample1)),
				ylab=deparse(substitute(sample2)) )
			chi2( sample1, sample2, bin=5 )
			kstwo( sample1, sample2 )

spectral.s		est.spectral.matrix( x1, y1, x2, y2, K=10 )
			spec.matrix( drifter1, drifter2, K=10 )
			plot.spec( f )
			sqrt.matrix( mat, inverse=F )
			plot.coh( f )
			sim.dep( f, N=500 )
			trace.stat( f )
			det.stat( f )
			indep.spec.test( f, K=10, method="determinant" )
			dep.spec.test( f, g, method="determinant" )
			wishart.sample( N, method="determinant" )

arima.s			ma1.sim( N, a1=diag(4) )
			ma2.sim( N, a1=diag(4), a2=diag(4) )
			ma.power( N=50, K=10, method="determinant" )
			
##############################################################################

make.vel.frame()
	transforms the "drifter" list into a big data frame. 
	The data are arranged in the increasing order of the numbers
	in the filenames.  Make sure to "read.vel" before making the
	data frame.   

read.vel() 
	reads real drifter data files into data frames and combines all 19
	data frames into a list called "drifter".  The function assigns
	individual data frames and the list to the .Data directory.
	
	There are six fields associated with each data file (or each
	data frame in the Splus data structure).  They are days, longitude,
	latitude, temperature, v1, and v2.  

	Note:  The real drifter data files are in the directory
	~awang/drifter/vel_data.  When read.vel() is called, the data
	files have to be moved to the directory where Splus is
	invoked.  

scan.drifter( filename )
	reads the simulated drifter files.  The output is a list of
	drifters, like "drifter" in the real data case.  Each drifter is
	arranged into a data frame with two fields, latitude and longitude.
	This function handles both the regular output file and the file with
	insertion from the simulation program.  

plot.drifter( drifter.list, scale=F )
	plots a list of drifters, like the output of read.vel() or
	scan.drifter().  If the option "scale" is TRUE, all the drifters
	will be plotted with the same scale.  

plot.traj( ..., inter=F )
	plots the trajectories of the given drifters and superimposes the
	parts of the drifters that intersect in time if the option "inter"
	is TRUE.

drifter.world()
	plots part of the world that all the drifters travel, which
	include the South Atlantic and Indian Oceans.

motion( starting=2446200, n=100, plot.all=F )
	plots all the existing drifters between the time frame starting from
	"starting for "n" days.  If the option plot.all is TRUE, all the
	drifters will be plotted.

wishart.ratio( N=100, r=2, K=10, method="determinant" )
	computes N trace statistics or determinant statistics of (W1 *
	W2^(-1)), where W1 and W2 have a Wishart distribution of dimension r
	and degrees of freedom K.

det( x )
	computes the determinant of the square matrix x.

myqq( sample1, sample2, xlab=deparse(substitute(sample1)),
			ylab=deparse(substitute(sample2)) )
	compares the distributions of two samples by qqplot.

chi2( sample1, sample2, bin=5 )
	performs chi-square two-sample test, but it does NOT work!  ;^b

kstwo( sample1, sample2 )
	performs Kolmogorov-Sirnov two-sample test and returns the
	p-value.

est.spectral.matrix( x1, y1, x2, y2, K=10 )
	computes the 4 by 4 spectral density matrix estimate for the
	given vector-valued series.

spec.matrix( drifter1, drifter2, K=10 )
	computes the 4 by 4 spectral density matrix estimate for the given
	drifters.  The drifters can be either two real drifters or tow
	dimulated drifters.

plot.spec( f )
	plots the modulus of the 4 by 4 spectral density matrix estimate at
	all calculated frequencies with the same scale in each panel.

sqrt.matrix( mat, inverse=F )
	computes the (inverse) square root of a matrix.

plot.coh( f )
	plots the coherency computed from the spectral density estimate at all
	calculated frequencies.

sim.dep( f, N=500 )
	simulates N by 4 data with the given spectral density matrix f.

trace.stat( f )
	computes the "trace" statistic of the spectral density matrix f.

det.stat( f )
	computes the "determinant" statistic of the spectral density matrix f.

indep.spec.test( f, K=10, method="determinant" )
	tests if either the trace or the determinant sample computed from f
	has the same distribution as the corresponding sample generated from
	Monte Carlo simulation.  It returns the p-value of the
	Kolmogorov-Sirmov two-sample test.

dep.spec.test( f, g, method="determinant" )
	tests if either the trace or the determinant sample computed from f
	has the same distribution as the corresponding sample computed from g.
	It returns the p-value of the Kolmogorov-Sirmov two-sample test.

wishart.sample <- function( N, method="determinant" )
	computes N trace statistics or determinant statistics of (W1 *
	W2^(-1)) from independent normal sample with K ranges from 2
	to 10.  

ma1.sim( N, a1=diag(4) )
	generates an N by 4 MA(1) series with coefficient matrix a1.  
	X(t) = e(t) + a1 * e(t-1), t = 1...N.

ma2.sim( N, a1=diag(4), a2=diag(4) )
	generates an N by 4 MA(2) series with coefficient matrices a1 and a2.
	X(t) = e(t) + a1 * e(t-1) + a2 * e(t-2), t = 1...N.

ma.power( N=50, K=10, np=2020 )
	computes the Prob{ accept H_0 | H_1 or H_0 } * N, i.e., the
	number of times accepting H_0 out of N trials, where H_0 : f_11 =
	f_22, and H_1 : f_11 != f_22.  The test statistics are the trace and
	the determinant of the matrix (W1 * (W2)^(-1)).  The independent test
	uses the histograms generated by Monte Carlo simulations.  They are
	wishart.trace and wishart.det.  The dependent test uses the histogram
	generated according to the estimate of the spectral density matrix.
