DEFINITION MODULE CycloStat;
(*******************************************************************
Module CycloStat (Version 1.0)
Copyright (c) 1995-2006 by Dimitrios Gyalistras, Andreas Fischlin
and ETH Zurich.
Purpose Estimate the parameters of cyclo-stationary,
multivariate AR(1)-type time series models.
Includes a procedure for Fourier analysis.
Remarks Missing data are recognized and returned as
UndefLONGREAL().
Programming
o Design
Dimitrios Gyalistras 13/07/1995
o Implementation
Dimitrios Gyalistras 13/07/1995
Andreas Fischlin 24/05/2002
ETH Zurich
Systems Ecology
CHN E 35.1
Universitaetstrasse 16
8092 Zurich
SWITZERLAND
URLs:
<mailto:RAMSES@env.ethz.ch>
<http://www.sysecol.ethz.ch>
<http://www.sysecol.ethz.ch/SimSoftware/RAMSES>
Last revision of definition: 24/05/2002 AF
*******************************************************************)
IMPORT Errors;
FROM LgMatrices IMPORT LMatrix, LVector;
(******************************************************************)
CONST (* result codes returned *)
allOk = Errors.allOk;
notDone = Errors.onlyAnInsert;
(*
Note: the dimensions of all LVector and LMatrix objects used below
are determined automatically via procedures LgMatrices.NElems,
LgMatrices.NRows and LgMatrices.NCols.
*)
CONST
MaxPeriodLgth = 512;
TYPE
LMatrixArr = ARRAY [1..MaxPeriodLgth] OF LMatrix;
(*--------------------------------------------------------------------------*)
(* Auxiliary procedures *)
(*--------------------------------------------------------------------------*)
PROCEDURE NBasisFuncs( periodLength, nHarmonics: INTEGER ): INTEGER;
(*
The number of basis functions used if "nHarmonics" pairs of
harmonics are considered for the expansion, and the period
length equals "periodLength".
*)
PROCEDURE EvalBasisFunc( bFuncNr, periodLength, phase: INTEGER ): LONGREAL;
(*
The value of the "bFuncNr"-th basis function at phase "phase"
given a period length "periodLength".
*)
(*--------------------------------------------------------------------------*)
(* Fourier Analysis *)
(*--------------------------------------------------------------------------*)
PROCEDURE Fourier( data : LVector; (* in: time series to analyze *)
periodLength : INTEGER; (* in: period length in time steps *)
nHarmonics : INTEGER; (* in: number of Sin/Cos pairs to evaluate *)
VAR mu : LONGREAL; (* out: long-term mean of ts *)
VAR A : LVector; (* out: amplitudes of Cosines *)
VAR B : LVector; (* out: amplitudes of Sines *)
VAR phase : LVector; (* out: phase [Sqrt(Ai^2+Bi^2)*Cos(omega*t-phase)] *)
VAR prdgr : LVector; (* out: periodogram *)
VAR vxpl : LVector; (* out: % variance expl. by freq. i *)
VAR resCode : INTEGER;
VAR errTxt : ARRAY OF CHAR );
PROCEDURE ReconstructCycle(
mu : LONGREAL; (* in: long-term mean of time series *)
A : LVector; (* in: amplitudes of Cosines *)
B : LVector; (* in: amplitudes of Sines *)
periodLength : INTEGER; (* in: period length in time steps *)
nHarmonics : INTEGER; (* in: number of Sin/Cos pairs *)
cycle : LVector; (* out: reconstructed cycle *)
VAR resCode : INTEGER;
VAR errTxt : ARRAY OF CHAR );
(*--------------------------------------------------------------------------*)
(* Fitting of cyclo-stationary models *)
(*--------------------------------------------------------------------------*)
(*
The following procedures estimate from a multivariate data set which
consists of "nTimeSeries" variables the "nTimeSeries x nTimeSeries"-
sized matrices "A(p)" and "B(p)" of the cyclo- stationary AR(1)-process
X(k+1) = A(p)*X(k) + B(p)*E(k)
Here "k" is the time step, "p = (k+periodLength) MOD periodLength"
is the phase within the assumed periodicity, "X" is the state
vector with E[X(i)] = 0 for all system states i, and "E" is a vector
of independent random components from a "nTimeSeries"-dimensional
normal distribution N(0,1).
The process matrices can be determined either separately for
each phase ("fixed-phase" model), or they can be determined
by first projecting the system states X onto the sub-space
spanned by the orthogonal basis functions
g1(k) = 1,
g2(k) = Cos( 2*PI*1*fo*k ), g3(k) = Sin( 2*PI*1*fo*k ),
g4(k) = Cos( 2*PI*2*fo*k ), g5(k) = Sin( 2*PI*2*fo*k ),
...
gM(k) = Cos( 2*PI*fmax*k ), gM+1(k) = Sin( 2*PI*fmax*k ),
where fo=(1/periodLength) and fmax<=fNyquist=periodLength/2.
("phase-averaged" model).
For fmax
PROCEDURE GetFixedPhaseCovs( data : LMatrix; (* in: rows=variables, columns=timepoints *)
periodLength : INTEGER; (* in: deterministic period in time steps *)
standardize : BOOLEAN; (* in: if true, cov-matrices are determined for the
time series scaled to unit standard deviation *)
VAR M : LMatrix; (* out: means for each phase, rows=phase *)
VAR S : LMatrix; (* out: standard deviations for each phase, rows=phase *)
VAR C0 : LMatrixArr; (* out: lag-0 covariances for each phase *)
VAR C1 : LMatrixArr; (* out: lag-1 covariances for each phase *)
VAR resCode : INTEGER;
VAR errTxt : ARRAY OF CHAR );
PROCEDURE FitFixedPhaseModel( periodLength : INTEGER;
C0 : LMatrixArr; (* in: lag-0 covariances for each phase *)
C1 : LMatrixArr; (* in: lag-1 covariances for each phase *)
A : LMatrixArr; (* out: A-matrices for each phase *)
B : LMatrixArr; (* out: B-matrices for each phase *)
VAR resCode : INTEGER;
VAR errTxt : ARRAY OF CHAR );
PROCEDURE FitPhaseAveragedModel( data : LMatrix; (* in: rows=variables, columns=timepoints *)
periodLength : INTEGER; (* in: deterministic period in time steps *)
standardize : BOOLEAN; (* in: if true, model is determined for
time series scaled to unit standard deviation *)
nHarmonics : INTEGER; (* in: number of harmonics to use, i.e. number of
Sin/Cos-pairs in addition to basis function "g1" *)
M : LMatrix; (* out: means for each phase *)
S : LMatrix; (* out: standard deviations for each phase *)
A : LMatrixArr; (* out: system matrices for each phase *)
B : LMatrixArr; (* out: error matrices for each phase *)
VAR badBEigVal : BOOLEAN; (* out: true if bad eigenvalue for B-matrix occurred *)
VAR resCode : INTEGER;
VAR errTxt : ARRAY OF CHAR );
END CycloStat.
RAMSES
Modules:
A-Z
Function
Layer
QuickRefs:
DM
AuxLib
AuxLibE
SciLib
EasyMW
MW
ISIS
RMSLib