|
|
|
|
||
|
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 fmaxPROCEDURE 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.
|
||
|
|
|