Running forecasts

In order to fit the simulation model to observations and generate forecasts, we need to:

  1. Choose prior distributions for \(x(t)\), \(y(t)\), and \(z(t)\);

  2. Define an input file for each observation model; and

  3. Record summary statistics such as predictive credible intervals for each observation model.

These changes are indicated by the highlighted lines in the following scenario definition:

An example scenario for generating forecasts for the Lorenz-63 system.
 1[components]
 2model = "pypfilt.examples.lorenz.Lorenz63"
 3time = "pypfilt.Scalar"
 4sampler = "pypfilt.sampler.LatinHypercube"
 5summary = "pypfilt.summary.HDF5"
 6
 7[time]
 8start = 0.0
 9until = 25.0
10steps_per_unit = 10
11summaries_per_unit = 10
12
13[prior]
14sigma = { name = "constant", args.value = 10 }
15rho = { name = "constant", args.value = 28 }
16beta = { name = "constant", args.value = 2.66667 }
17x = { name = "uniform", args.loc = -5, args.scale = 10 }
18y = { name = "uniform", args.loc = -5, args.scale = 10 }
19z = { name = "uniform", args.loc = -5, args.scale = 10 }
20
21[observations.x]
22model = "pypfilt.examples.lorenz.ObsLorenz63"
23file = "lorenz63-x.ssv"
24
25[observations.y]
26model = "pypfilt.examples.lorenz.ObsLorenz63"
27file = "lorenz63-y.ssv"
28
29[observations.z]
30model = "pypfilt.examples.lorenz.ObsLorenz63"
31file = "lorenz63-z.ssv"
32
33[summary.tables]
34forecasts.component = "pypfilt.summary.PredictiveCIs"
35forecasts.credible_intervals = [50, 60, 70, 80, 90, 95]
36
37[filter]
38particles = 500
39prng_seed = 2001
40history_window = -1
41resample.threshold = 0.25
42
43[scenario.forecast]

Observations will be read from these files when generating forecasts with pypfilt.forecast():

def run_lorenz63_forecast():
    scenario_file = 'lorenz63_forecast.toml'
    with open(scenario_file, 'w') as f:
        f.write(pypfilt.examples.lorenz.lorenz63_forecast_toml())
    instances = list(pypfilt.load_instances(scenario_file))
    instance = instances[0]

    # Run a forecast from t = 20.
    forecast_time = 20
    context = instance.build_context()
    return pypfilt.forecast(context, [forecast_time], filename=None)

If you pass a filename to pypfilt.forecast(), all of the summary tables for the estimation pass and each forecasting pass will be saved to that file, as HDF5 data sets.

Note

HDF5 is a file format that allows you to store lots of data tables and related metadata in a single file, and to load these data tables as if they were NumPy arrays. All of the summary tables recorded by pypfilt are NumPy structured arrays. You can explore HDF5 files with the h5py package, which makes it easy to load and store data tables.