cats_v312.tex

最大似然估计算法

\documentclass[12pt]{article}
\usepackage{times}

\begin{document}

\title{Create and Analyse Time Series :  CATS software V3.1.1}
\author{Simon Williams}
\maketitle

\section{Introduction}

{\bf cats} is a program to use Maximum Likelihood Estimation to fit a multi-parameter
model to a time series (such as continuous GPS).  The routine solves
for all parameters simultaneously, but in two parts to increase speed;
the linear part includes an offset and slope,  the possibility of
abrupt steps (as from earthquakes) and sinusoidal terms (for example annual and
semiannual terms), while the non-linear part solves for
several specific noise models and combinations.

\begin{itemize}

\item
White noise.

\item 
Power-law noise (otherwise known as fractionally integrated noise)

\item 
First-Order Gauss Markov noise (equivalent to autoregressive AR(1) noise)

\item
Band-pass noise (See {\it Langbein} [2004])

\item
Generalized Gauss Markov noise (See {\it Langbein} [2004])

\item
Variable white noise (i.e. using the formal errors)

\item
Step-variable white noise (i.e. a change in the scale of white noise between two epochs)

\item
Time-variable white noise (i.e. an exponential decay in the scale of the white noise)

\end{itemize}

The background to this program can be found in the following papers, 
{\it Langbein and Johnson} [1997], {\it Zhang et al.} [1997], 
{\it Mao et al.} [1999], {\it Williams} [2003], {\it Williams et al.} [2004] and {\it Langbein} [2004]. 
The program is command-line orientated, so it does not require any other files
except for the time series file in order to run. 
For the power-law noise models, {\bf cats} uses the fractional difference
method [{\it Hosking}, 1981] for creating the covariance matrix. 
{\it Zhang et al.} [1997] and {\it Mao et al.} [1999] used an approximate covariance
matrix for flicker noise derived from an algorithm described by {\it Gardner} [1978]
for simulating such noise (H. Johnson, pers. comm., 1996). The constants in this matrix 
were chosen so that the power spectrum of flicker noise and
random walk noise cross at a period of one year when both have a sampling interval of one day
and equal amplitude ($b_{-1} = b_{-2}$). This earlier covariance matrix is not exactly the same as that
derived from the above transformation matrix. The major difference is in the scaling of the amplitude.
The scaling between the ``new'' and ``old'' amplitudes can be estimated from their respective power
spectrum equations. That is
\begin{eqnarray}
P_{old} & = & \frac{b_{old}^2 f^{-1}}{2 \pi^2} \\ \nonumber
P_{new} & = & \frac{b_{new}^2 f^{-1}}{\pi \sqrt{f_s \times 24 \times 60
\times 60 \times 365.249}}
\end{eqnarray}
therefore
\begin{equation}
b_{new} = \frac{ (f_s \times 24 \times 60 \times 60 \times 365.249)^{1/4} }
{\sqrt{2 \pi}} \; b_{old}.
\end{equation}
When the sampling frequency is once per day then
\begin{equation}
b_{new} = 1.7440 \; b_{old}.
\end{equation}
Therefore, values quoted for the amplitude of flicker noise in papers that use the {\it Zhang et al.} [1997] 
covariance matrix must be scaled before being used in any studies using the new matrix. 

In the formulisation used in this program, the system is scaled so
that for any spectral index, the power spectra cross over at the same
frequency given the same sampling frequency and the same size $\sigma$. 
For daily sampling, the cross over frequency is about 120 days. This also
means that for spectral indices other than random walk and white noise
the scaling parameter changes with different sampling frequencies. This
must be taken into account when comparing results from daily solutions
to those from weekly solutions or if the data has been decimalised. For instance,
given flicker noise and a sampling interval of 2 days, the flicker noise
parameter will be $2^{1/4}$ times smaller than the daily value. John Langbein
uses a different scheme again so that flicker noise amplitudes from
this program will be 4.3717 times larger than his results. 

\section{The Time Series File}

There are currently two accepted formats for entry; a ``cats'' file
and a ``psmsl'' file. In the future, other formats may become available. Users
are free to suggest new formats (writing the code would be even better -
see read$\_$series.c in the lib directory for examples)

\subsection{``cats'' file description}

The data file is fairly free format consisting of two parts, header
information and the time series. The header information is basically 
a list of parameters relevant to the GPS site the time series comes
from. The header information (plus any line you want to eliminate
from the processing) starts with a $\#$ symbol. The only header information
relevant to {\bf cats} is the list of offsets. The data set consists
of seven columns (no particular format) corresponding to, in order,
time (in decimal years) north, east, up, north error, east error, up error.
The positions are, by default, assumed to be in metres, however this
can be altered with the -Sscale$\_$factor option. At present the position
uncertainties are not used in the program. In the future the individual formal errors
may be used to define a variable white noise component. However this slows the algorithm
down considerably when compared to assuming the white noise is the
same at each epoch. If, for instance, there was only one component of
information required for processing (for example tide gauge data) then the other columns
must still be present but the program can be told not to process those
columns by using the -\--columns (-C) option described below. An example of a file
following this format is shown below.

{\scriptsize
\begin{verbatim}
# Site : vyas
# X :  -2483507.1216
# Y :  -4672361.9853
# Z :   3549320.5426
# Latitude : 34.030915
# Longitude : -117.992047
# Height : 7.856557
# offset 1999.79178082 7
1998.5736   -0.02370    0.04070    0.00290  0.0009  0.0009  0.0037
1998.5763   -0.02360    0.04190    0.00410  0.0009  0.0009  0.0030
1998.5791   -0.02180    0.04140    0.00650  0.0010  0.0009  0.0038
1998.5818   -0.02270    0.04100   -0.00130  0.0009  0.0009  0.0038
1998.5845   -0.02330    0.04070    0.00050  0.0009  0.0009  0.0037
1998.5873   -0.02300    0.04050   -0.00020  0.0009  0.0009  0.0038
1998.5900   -0.02320    0.04090    0.00560  0.0009  0.0009  0.0037
1998.5928   -0.02390    0.04030   -0.00180  0.0010  0.0009  0.0039
1998.5955   -0.02370    0.04020    0.00230  0.0009  0.0009  0.0037
1998.5983   -0.02320    0.04150    0.00480  0.0010  0.0009  0.0038
1998.6010   -0.02390    0.04020   -0.00620  0.0010  0.0009  0.0039
1998.6036   -0.02240    0.03970   -0.00140  0.0010  0.0009  0.0038
1998.6091   -0.02190    0.03990   -0.00190  0.0010  0.0009  0.0030
\end{verbatim}
}

If there any known discontinuities or offsets in the time series then
these should be specified in the header of the file in the following
format
{\scriptsize
\begin{verbatim}
# offset decimal_date component_code
\end{verbatim}
}
The component code is calculated in a manner similar to file permissions
on a unix-style system. The component code is the sum of the components
in which the offsets appears with the values north (4), east (2) and
vertical (1). This gives an exclusive number between 0 and 7 for all
combinations of affected components. For example for an offset occurring
in the north and up components the component code is 5. For all components
the value is 7. The offset line is repeated for as many offsets as there
are thought to be in the time series. The offsets need not be in chronological
order nor do you even have to ensure that they are within the time span of
the data set. The program will sort the offsets correctly when processing the
file.

\subsection{``psmsl'' file description}

This type of data file is again very simplistic and is described on the Permanent
Service for Mean Sea Level web site (http://www.pol.ac.uk/psmsl/datainfo). It
basically consists of one two main columns. Column one is the time stamp and
column two is the data. {\bf Cats}, since it is written in C, is not so specific
about the exact format, it just tries to read 3 columns on each line. If there
are two columns (time and data) then it uses that epoch. An example of a data file
is given below.
{\scriptsize
\begin{verbatim}
   1964.625  7036
   1964.708  7103
   1964.792  7057
   1964.875  7139
   1964.958  7091
   1965.042  7072
   1965.125  6834
   1965.208  6935
   1965.292  6941
   1965.375  6920
   1965.458  7011
   1965.542  6974 12
   1965.625  7026  8
   1965.708  7078  3
   1965.792  7078  3
   1965.875  7075
   1965.958  7115
   1966.042  6990
   1966.125  7087
   1966.208  7005
\end{verbatim}
}

As for the ``cats'' format, this file format has been adapted so that it can read in offsets and
ignore lines that begin with the $\#$ symbol.

\section{Command Line Options}

The command line options now come in two ``flavours'' depending on whether your machine
architecture can handle long options or not. In the following section the long options
will be described with the short option listed in brackets. There is a slight difference
between the version that can accept long options to the version that cannot. The long option
version would look something like this
{\scriptsize
\begin{verbatim}
cats vyas.neu --sinusoid 1y1 --verbose --columns 4 --output vyas_all.mle
\end{verbatim}
}
or
{\scriptsize
\begin{verbatim}
cats vyas.neu --sinusoid=1y1 --verbose --columns=4 --output=vyas_all.mle
\end{verbatim}
}
or
{\scriptsize
\begin{verbatim}
cats vyas.neu -A 1y1 -V -C4 -Ovyas_all.mle
\end{verbatim}
}
It can have a mixture of long and short. If using long options, then options with additional
parameters should either have a space or an equal sign between the option and the parameter/s.
If you choose to use short options then the additional parameter should follow on from the
option (e.g. -C4) or there should be a space between (e.g. -C 4).
For the version that only accepts short options then the options should always take the form
{\scriptsize
\begin{verbatim}
cats vyas.neu -A1y1 -V -C4 -Ovyas_all.mle
\end{verbatim}
}
with no gaps between option and the additional parameters.

\begin{description}

\item[-\--model (-M)] A description of the stochastic model. The first two characters
in the string indicate the type of stochastic model to use. These include {\bf pl} for power-law noise,
{\bf gm} for first-order gauss markov, {\bf wh} for white noise, {\bf bp} for band-pass noise, {\bf vw} for
variable white noise, {\bf gg} for generalised gauss-markov, and {\bf sw} for step-variable white noise.

For the power-law noise model you can fix the model to a specific spectral index by appending an additional
parameter as follows 
\begin{verbatim}
--model pl:k-1
\end{verbatim}
The above example would fix the spectral index to -1 (flicker noise). If a spectral index is not
specified then the program will also try to solve for the spectral index.

For the first-order gauss markov model you can fix the model to a specific $\beta$ (roughly
equivalent to the cross-over frequency) by appending an additional parameter as follows
\begin{verbatim}
--model gm:b23.0
\end{verbatim}
The above example would fix $\beta$ to 23.0. If $\beta$ is not 
specified then the program will also try to solve for $\beta$.

For the band-pass model there are three parameters $f_{central}$ (c), width (w) 
and the number of poles (p). You can fix any of these parameters as follows
\begin{verbatim}
--model bp:c1w1p1
\end{verbatim}
The band-pass model is defined slightly differently to that in {\it Langbein} [2004].
Instead of the parameters $f_l$ and $f_h$ to define the limits of the passband, we
use $f_{c}$, the central frequency, and $width$, the width of the passband. The two
sets of parameters are related using the following
\begin{eqnarray}
f_l & = & \frac{f_c}{width} \\ \nonumber
f_h & = & f_c  width \\ \nonumber
\end{eqnarray}
or alternatively
\begin{eqnarray}
f_c & = & \sqrt{f_h f_l} \\ \nonumber
width & = & \sqrt{\frac{f_l}{f_h}} \\ \nonumber
\end{eqnarray}

The variable white noise model uses the formal errors recorded in the data file (only
for the {\bf cats} files) and solves for a scale parameter (as opposed to a noise
amplitude for white noise). The are no other parameters for this model.

The generalised-gauss-markov model is described in {\it Langbein} [2004]. There are
two parameters for this model, equivalent to the spectral index of the power-law noise