You can run this notebook in a live session or view it on Github.

Calculate skill of a MJO Index of SubX model GEOS_V2p1 as function of weekly lead time

In this example, we demonstrate:

1. How to remotely access data from the Subseasonal Experiment (SubX) hindcast database and set it up to be used in climpred.

2. How to calculate the Anomaly Correlation Coefficient (ACC) using weekly data with climpred

3. How to calculate and plot historical forecast skill of the real-time multivariate MJO (RMM) indices as function of lead time.

The Subseasonal Experiment (SubX)

Further information on SubX is available from Pegion et al. 2019 and the SubX project website.

The SubX public database is hosted on the International Research Institute for Climate and Society (IRI) data server http://iridl.ldeo.columbia.edu/SOURCES/.Models/.SubX/.

Since the SubX data server is accessed via this notebook, the time for the notebook to run may is several minutes and will vary depending on the speed that data can be downloaded. This is a large dataset, so please be patient. If you prefer to download SubX data locally, scripts are available from kpegion/SubX.

Definitions

RMM

Two indices (RMM1 and RMM2) are used to represent the MJO. Together they define the MJO based on 8 phases and can represent both the phase and amplitude of the MJO [Wheeler and Hendon, 2004]. This example uses the observed RMM1 provided by Matthew Wheeler at the Center for Australian Weather and Climate Research. It is the version of the indices in which interannual variability has not been removed.

Skill of RMM

Traditionally, the skill of the RMM is calculated as a bivariate correlation encompassing the skill of the two indices together [Gottschalck et al., 2010, Rashid et al., 2011]. Currently, climpred does not have the functionality to calculate the bivariate correlation, thus the anomaly correlation coefficient for RMM1 index is calculated here as a demonstration. The bivariate correlation metric will be added in a future version of climpred.

# linting
%load_ext nb_black
%load_ext lab_black

import warnings

import matplotlib.pyplot as plt

import xarray as xr
import pandas as pd
import numpy as np

from climpred import HindcastEnsemble
import climpred

warnings.filterwarnings("ignore")

Read the observed RMM Indices

obsds = climpred.tutorial.load_dataset("RMM-INTERANN-OBS")[["rmm1"]]
obsds

<xarray.Dataset>
Dimensions:  (time: 15613)
Coordinates:
  * time     (time) datetime64[ns] 1974-06-03 1974-06-04 ... 2017-07-24
Data variables:
    rmm1     (time) float64 1.867 1.868 1.924 1.574 ... 0.8918 0.9023 1.028

xarray.Dataset

• Dimensions:
  • time: 15613
• Coordinates: (1)
  • time
    (time)
    datetime64[ns]
    1974-06-03 ... 2017-07-24

    array(['1974-06-03T00:00:00.000000000', '1974-06-04T00:00:00.000000000',
           '1974-06-05T00:00:00.000000000', ..., '2017-07-22T00:00:00.000000000',
           '2017-07-23T00:00:00.000000000', '2017-07-24T00:00:00.000000000'],
          dtype='datetime64[ns]')

• Data variables: (1)
  • rmm1
    (time)
    float64
    ...

    array([1.86732 , 1.86838 , 1.92386 , ..., 0.891819, 0.902347, 1.027518])

• Attributes: (0)

Read the SubX RMM1 data for the GMAO-GEOS_V2p1 model from the SubX data server. It is important to note that the SubX data contains weekly initialized forecasts where the init day varies by model. SubX data may have all NaNs for initial dates in which a model does not make a forecast, thus we apply dropna over the S=init dimension when how=all data for a given S=init is missing. This can be slow, but allows the rest of the calculations to go more quickly.

Note that we ran the dropna operation offline and then uploaded the post-processed SubX dataset to the climpred-data repo for the purposes of this demo. This is how you can do this manually:

url = 'http://iridl.ldeo.columbia.edu/SOURCES/.Models/.SubX/.GMAO/.GEOS_V2p1/.hindcast/.RMM/.RMM1/dods/'
fcstds = xr.open_dataset(url, decode_times=False, chunks={'S': 1, 'L': 45}).dropna(dim='S',how='all')

fcstds = climpred.tutorial.load_dataset("GMAO-GEOS-RMM1", decode_times=False)
fcstds

<xarray.Dataset>
Dimensions:  (S: 510, M: 4, L: 45)
Coordinates:
  * S        (S) float32 1.424e+04 1.425e+04 1.426e+04 ... 2.044e+04 2.045e+04
  * M        (M) float32 1.0 2.0 3.0 4.0
  * L        (L) float32 0.5 1.5 2.5 3.5 4.5 5.5 ... 40.5 41.5 42.5 43.5 44.5
Data variables:
    RMM1     (S, M, L) float32 -0.03237 0.1051 0.1946 ... -0.1499 0.08326 0.1389
Attributes:
    Conventions:  IRIDL

xarray.Dataset

• Dimensions:
  • S: 510
  • M: 4
  • L: 45
• Coordinates: (3)
  • S
    (S)
    float32
    1.424e+04 1.425e+04 ... 2.045e+04

    long_name :

    Start Time

    calendar :

    standard

    standard_name :

    forecast_reference_time

    pointwidth :

    0.0

    gridtype :

    0

    units :

    days since 1960-01-01

    array([14245., 14250., 14255., ..., 20439., 20444., 20449.], dtype=float32)

  • M
    (M)
    float32
    1.0 2.0 3.0 4.0

    long_name :

    Ensemble Member

    standard_name :

    realization

    pointwidth :

    1.0

    gridtype :

    0

    units :

    ids

    array([1., 2., 3., 4.], dtype=float32)

  • L
    (L)
    float32
    0.5 1.5 2.5 3.5 ... 42.5 43.5 44.5

    long_name :

    Lead

    standard_name :

    forecast_period

    pointwidth :

    1.0

    gridtype :

    0

    units :

    days

    array([ 0.5,  1.5,  2.5,  3.5,  4.5,  5.5,  6.5,  7.5,  8.5,  9.5, 10.5, 11.5,
           12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5, 19.5, 20.5, 21.5, 22.5, 23.5,
           24.5, 25.5, 26.5, 27.5, 28.5, 29.5, 30.5, 31.5, 32.5, 33.5, 34.5, 35.5,
           36.5, 37.5, 38.5, 39.5, 40.5, 41.5, 42.5, 43.5, 44.5], dtype=float32)

• Data variables: (1)
  • RMM1
    (S, M, L)
    float32
    ...

    pointwidth :

    0.0

    long_name :

    RMM1

    units :

    unitless

    array([[[-0.032366,  0.105126, ...,  1.778099,  1.803143],
            [ 0.006499,  0.158598, ...,  0.826527,  0.743905],
            [ 0.014169,  0.147751, ...,  1.170505,  1.230083],
            [-0.014136,  0.142899, ...,  1.239202,  1.26709 ]],
    
           [[ 0.222802,  0.70556 , ...,  2.06537 ,  1.935836],
            [ 0.192239,  0.610177, ...,  1.658914,  1.51989 ],
            [ 0.200765,  0.662475, ...,  1.159152,  1.093163],
            [ 0.203409,  0.631355, ..., -0.66206 , -0.785497]],
    
           ...,
    
           [[ 0.167977,  0.093171, ...,  0.455726,  0.463756],
            [ 0.300361,  0.193978, ..., -0.853684, -0.762082],
            [ 0.223478,  0.095381, ..., -0.784702, -0.75755 ],
            [ 0.228956,  0.144475, ..., -0.866805, -0.977645]],
    
           [[-0.648854, -1.236488, ...,  0.467468,  0.470097],
            [-0.566065, -1.12568 , ...,  0.261335,  0.253228],
            [-0.638935, -1.209603, ..., -0.003442, -0.231822],
            [-0.603829, -1.186343, ...,  0.083258,  0.138907]]], dtype=float32)

• Attributes: (1)

  Conventions :

  IRIDL

The SubX data dimensions correspond to the following climpred dimension definitions: X=lon,L=lead,Y=lat,M=member, S=init. We will rename the dimensions to their climpred names.

fcstds = fcstds.rename({"S": "init", "L": "lead", "M": "member", "RMM1": "rmm1"})

Let’s make sure that the lead dimension is set properly for climpred. SubX data stores leads as 0.5, 1.5, 2.5, etc, which correspond to 0, 1, 2, … days since initialization. We will change the lead to be integers starting with zero.

fcstds["lead"] = (fcstds["lead"] - 0.5).astype("int")

Now we need to make sure that the init dimension is set properly for climpred. We use an xarray convenience function to decode it into the proper CFTime calendar. It can detect that this is on a 360 day calendar.

fcstds = xr.decode_cf(fcstds, decode_times=True)

Make Weekly/7-day Averages

fcstweekly = fcstds.rolling(lead=7, center=False).mean().dropna(dim="lead")
obsweekly = obsds.rolling(time=7, center=False).mean().dropna(dim="time")
print(fcstweekly)
print()
print(obsweekly)

<xarray.Dataset>
Dimensions:  (init: 510, member: 4, lead: 39)
Coordinates:
  * init     (init) datetime64[ns] 1999-01-01 1999-01-06 ... 2015-12-27
  * member   (member) float32 1.0 2.0 3.0 4.0
  * lead     (lead) int64 6 7 8 9 10 11 12 13 14 ... 36 37 38 39 40 41 42 43 44
Data variables:
    rmm1     (init, member, lead) float32 0.1347 0.1095 ... 0.1404 0.07606
Attributes:
    Conventions:  IRIDL

<xarray.Dataset>
Dimensions:  (time: 15456)
Coordinates:
  * time     (time) datetime64[ns] 1974-06-09 1974-06-10 ... 2017-07-24
Data variables:
    rmm1     (time) float64 1.336 1.107 0.9046 0.695 ... 0.6265 0.673 0.7352

fcstweekly = fcstweekly.isel(lead=slice(None, None, 7)).assign_coords(
    lead=np.arange(1, 1 + fcstweekly.lead.size // 7 + 1)
)

climpred requires that lead dimension has an attribute called units indicating what time units the lead is assocated with. Options are: years, seasons, months, weeks, pentads, days, hours, minutes or seconds. The lead units are weeks.

fcstweekly["lead"].attrs = {"units": "weeks"}

Create the HindcastEnsemble object and HindcastEnsemble.add_observations().

hindcast = HindcastEnsemble(fcstweekly)
hindcast = hindcast.add_observations(obsweekly)
hindcast

climpred.HindcastEnsemble

<Initialized Ensemble>
Dimensions:     (init: 510, member: 4, lead: 6)
Coordinates:
  * init        (init) object 1999-01-01 00:00:00 ... 2015-12-27 00:00:00
  * member      (member) float32 1.0 2.0 3.0 4.0
  * lead        (lead) int64 1 2 3 4 5 6
    valid_time  (lead, init) object 1999-01-08 00:00:00 ... 2016-02-07 00:00:00
Data variables:
    rmm1        (init, member, lead) float32 0.1347 -0.7031 ... 0.2992 0.4478
Attributes:
    Conventions:  IRIDL

Initialized Ensemble

• Dimensions:
  • init: 510
  • member: 4
  • lead: 6
• Coordinates: (4)
  • init
    (init)
    object
    1999-01-01 00:00:00 ... 2015-12-...

    standard_name :

    forecast_reference_time

    long_name :

    Initialization

    description :

    The forecast reference time in NWP is the "data time", the time of the analysis from which the forecast was made. It is not the time for which the forecast is valid; the standard name of time should be used for that time.

    array([cftime.DatetimeProlepticGregorian(1999, 1, 1, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(1999, 1, 6, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(1999, 1, 11, 0, 0, 0, 0, has_year_zero=True),
           ...,
           cftime.DatetimeProlepticGregorian(2015, 12, 17, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2015, 12, 22, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2015, 12, 27, 0, 0, 0, 0, has_year_zero=True)],
          dtype=object)

  • member
    (member)
    float32
    1.0 2.0 3.0 4.0

    long_name :

    Member

    standard_name :

    realization

    pointwidth :

    1.0

    gridtype :

    0

    units :

    ids

    description :

    Realization is used to label a dimension that can be thought of as a statistical sample, e.g., labelling members of a model ensemble.

    array([1., 2., 3., 4.], dtype=float32)

  • lead
    (lead)
    int64
    1 2 3 4 5 6

    units :

    weeks

    standard_name :

    forecast_period

    long_name :

    Lead

    description :

    Forecast period is the time interval between the forecast reference time and the validity time. A period is an interval of time, or the time-period of an oscillation.

    array([1, 2, 3, 4, 5, 6])

  • valid_time
    (lead, init)
    object
    1999-01-08 00:00:00 ... 2016-02-...

    long_name :

    validity time

    standard_name :

    time

    description :

    time for which the forecast is valid

    calculate :

    init + lead

    amip :

    time

    array([[cftime.DatetimeProlepticGregorian(1999, 1, 8, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(1999, 1, 13, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(1999, 1, 18, 0, 0, 0, 0, has_year_zero=True),
            ...,
            cftime.DatetimeProlepticGregorian(2015, 12, 24, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(2015, 12, 29, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(2016, 1, 3, 0, 0, 0, 0, has_year_zero=True)],
           [cftime.DatetimeProlepticGregorian(1999, 1, 15, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(1999, 1, 20, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(1999, 1, 25, 0, 0, 0, 0, has_year_zero=True),
            ...,
            cftime.DatetimeProlepticGregorian(2015, 12, 31, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(2016, 1, 5, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(2016, 1, 10, 0, 0, 0, 0, has_year_zero=True)],
           [cftime.DatetimeProlepticGregorian(1999, 1, 22, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(1999, 1, 27, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(1999, 2, 1, 0, 0, 0, 0, has_year_zero=True),
            ...,
            cftime.DatetimeProlepticGregorian(2016, 1, 7, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(2016, 1, 12, 0, 0, 0, 0, has_year_zero=True),
    ...
            cftime.DatetimeProlepticGregorian(1999, 2, 8, 0, 0, 0, 0, has_year_zero=True),
            ...,
            cftime.DatetimeProlepticGregorian(2016, 1, 14, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(2016, 1, 19, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(2016, 1, 24, 0, 0, 0, 0, has_year_zero=True)],
           [cftime.DatetimeProlepticGregorian(1999, 2, 5, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(1999, 2, 10, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(1999, 2, 15, 0, 0, 0, 0, has_year_zero=True),
            ...,
            cftime.DatetimeProlepticGregorian(2016, 1, 21, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(2016, 1, 26, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(2016, 1, 31, 0, 0, 0, 0, has_year_zero=True)],
           [cftime.DatetimeProlepticGregorian(1999, 2, 12, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(1999, 2, 17, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(1999, 2, 22, 0, 0, 0, 0, has_year_zero=True),
            ...,
            cftime.DatetimeProlepticGregorian(2016, 1, 28, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(2016, 2, 2, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(2016, 2, 7, 0, 0, 0, 0, has_year_zero=True)]],
          dtype=object)

• Data variables: (1)
  • rmm1
    (init, member, lead)
    float32
    0.1347 -0.7031 ... 0.2992 0.4478

    pointwidth :

    0.0

    long_name :

    RMM1

    units :

    unitless

    array([[[ 0.13466184, -0.70308524, -0.6665866 ,  0.53768194,
              1.5258206 ,  1.3847538 ],
            [ 0.2424833 ,  0.04300196, -0.04560124,  1.5524023 ,
              1.9928336 ,  1.306517  ],
            [ 0.29170546, -0.40440923,  0.289992  ,  0.902887  ,
              1.1658235 ,  1.5833153 ],
            [ 0.19159952, -1.0362829 , -0.41041547,  0.9215896 ,
              1.9797173 ,  1.7849027 ]],
    
           [[ 0.13745967, -0.66343164,  0.7237875 ,  0.8539322 ,
              1.4124439 ,  2.0646007 ],
            [ 0.05160916, -1.2358718 ,  0.07569085,  1.174731  ,
              0.74607867,  1.4797692 ],
            [ 0.15610577, -1.2223252 ,  0.42348805,  0.9170564 ,
              1.0033319 ,  1.0180309 ],
            [ 0.06007081, -0.943257  , -0.02319996,  1.3574636 ,
              2.5858319 ,  0.9067065 ]],
    
           [[-0.98613304,  0.5226586 ,  1.2636639 ,  0.92440665,
              1.3922706 ,  0.17960559],
    ...
            [ 0.5657733 , -1.343183  , -3.0290492 , -2.6482596 ,
             -1.8979372 , -0.8949425 ]],
    
           [[-0.07046235, -2.089687  , -3.5156634 , -2.9199138 ,
             -2.470964  , -0.8262414 ],
            [ 0.0170458 , -1.6903168 , -2.8156073 , -2.3722534 ,
             -0.88731194, -0.6438968 ],
            [-0.01062303, -2.2151313 , -3.0623872 , -3.8109653 ,
             -2.7191453 , -2.234973  ],
            [ 0.00703398, -1.7830768 , -2.8474169 , -2.4715648 ,
             -2.3030205 , -1.6926976 ]],
    
           [[-1.9386648 , -2.7401743 , -2.6152925 , -2.758163  ,
             -1.4414849 ,  0.1724483 ],
            [-1.7458452 , -3.2379034 , -2.7567887 , -1.8445466 ,
             -1.5084535 ,  0.1758536 ],
            [-1.8588039 , -3.2608886 , -3.4155226 , -1.1434904 ,
             -0.96848327,  0.2189652 ],
            [-1.8091087 , -2.9087188 , -3.6594112 , -3.0632272 ,
              0.29924652,  0.44781074]]], dtype=float32)

• Attributes: (1)

  Conventions :

  IRIDL

<Observations>
Dimensions:  (time: 15456)
Coordinates:
  * time     (time) object 1974-06-09 00:00:00 ... 2017-07-24 00:00:00
Data variables:
    rmm1     (time) float64 1.336 1.107 0.9046 0.695 ... 0.6265 0.673 0.7352

Observations

• Dimensions:
  • time: 15456
• Coordinates: (1)
  • time
    (time)
    object
    1974-06-09 00:00:00 ... 2017-07-...

    long_name :

    time

    standard_name :

    time

    array([cftime.DatetimeProlepticGregorian(1974, 6, 9, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(1974, 6, 10, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(1974, 6, 11, 0, 0, 0, 0, has_year_zero=True),
           ...,
           cftime.DatetimeProlepticGregorian(2017, 7, 22, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2017, 7, 23, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2017, 7, 24, 0, 0, 0, 0, has_year_zero=True)],
          dtype=object)

• Data variables: (1)
  • rmm1
    (time)
    float64
    1.336 1.107 0.9046 ... 0.673 0.7352

    array([1.33585898, 1.10670328, 0.90462442, ..., 0.62653254, 0.6730053 ,
           0.73522198])

• Attributes: (0)

Calculate the Anomaly Correlation Coefficient (ACC) climpred.metrics._pearson_r().

skill = hindcast.verify(
    metric="acc", comparison="e2o", alignment="maximize", dim="init"
)
skill

<xarray.Dataset>
Dimensions:  (lead: 6)
Coordinates:
  * lead     (lead) int64 1 2 3 4 5 6
    skill    <U11 'initialized'
Data variables:
    rmm1     (lead) float64 0.948 0.8719 0.7265 0.5658 0.4504 0.3607
Attributes:
    Conventions:                   IRIDL
    prediction_skill_software:     climpred https://climpred.readthedocs.io/
    skill_calculated_by_function:  HindcastEnsemble.verify()
    number_of_initializations:     510
    number_of_members:             4
    alignment:                     maximize
    metric:                        pearson_r
    comparison:                    e2o
    dim:                           init
    reference:                     []

xarray.Dataset

• Dimensions:
  • lead: 6
• Coordinates: (2)
  • lead
    (lead)
    int64
    1 2 3 4 5 6

    units :

    weeks

    standard_name :

    forecast_period

    long_name :

    Lead

    description :

    Forecast period is the time interval between the forecast reference time and the validity time. A period is an interval of time, or the time-period of an oscillation.

    array([1, 2, 3, 4, 5, 6])

  • skill
    ()
    <U11
    'initialized'

    description :

    new dimension prediction skill of initialized and reference forecasts created by .verify() or .bootstrap()

    documentation :

    https://climpred.readthedocs.io/en/v2.1.7.dev40+ga522c2d.d20211213/reference_forecast.html

    array('initialized', dtype='<U11')

• Data variables: (1)
  • rmm1
    (lead)
    float64
    0.948 0.8719 ... 0.4504 0.3607

    pointwidth :

    0.0

    long_name :

    RMM1

    units :

    None

    array([0.94796002, 0.87192995, 0.72647792, 0.56582236, 0.45044375,
           0.36065706])

• Attributes: (10)

  Conventions :

  IRIDL

  prediction_skill_software :

  climpred https://climpred.readthedocs.io/

  skill_calculated_by_function :

  HindcastEnsemble.verify()

  number_of_initializations :

  510

  number_of_members :

  4

  alignment :

  maximize

  metric :

  pearson_r

  comparison :

  e2o

  dim :

  init

  reference :

  []

Plot the skill as a function of lead time

skill.rmm1.to_series().plot(kind="bar")
plt.title("GMAO_GOES_V2p1 RMM1 Skill")
plt.xlabel("Lead Time (Weeks)")
plt.ylabel("ACC")
plt.ylim(0.0, 1.0)
plt.show()

References

[1]

J. Gottschalck, M. Wheeler, K. Weickmann, F. Vitart, N. Savage, H. Lin, H. Hendon, D. Waliser, K. Sperber, M. Nakagawa, C. Prestrelo, M. Flatau, and W. Higgins. A Framework for Assessing Operational Madden–Julian Oscillation Forecasts: A CLIVAR MJO Working Group Project. Bulletin of the American Meteorological Society, 91(9):1247–1258, September 2010. doi:10/fcftn7.

[2]

Harun A. Rashid, Harry H. Hendon, Matthew C. Wheeler, and Oscar Alves. Prediction of the Madden–Julian oscillation with the POAMA dynamical prediction system. Climate Dynamics, 36(3):649–661, February 2011. doi:10/b992k4.

[3]

Matthew C. Wheeler and Harry H. Hendon. An All-Season Real-Time Multivariate MJO Index: Development of an Index for Monitoring and Prediction. Monthly Weather Review, 132(8):1917–1932, August 2004. doi:10/d489vd.