ETHZ_Logo RAMSES_Logo_Right   RAMSES   RAMSES_Logo_Left Systems Ecology  
Start    search button      Modules:   A-Z   Function   Layer        QuickRefs:   DM   AuxLib   AuxLibE   SciLib   EasyMW   MW   ISIS   RMSLib

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.

  Contact RAMSES@env.ethz.ch Last updated: 25-Jul-2011 [Top of page]