Chapter 8.3 Forecast Experiment and Predictive Skill

Learning Objectives

After completing this chapter, you should be able to

initialize a forecast using observations,
perform forecast and hindcast experiments in a perfect-model framework,
compare forecasts with the true climate state,
evaluate forecast quality using RMSE, ACC, and RMSE Skill Score (RMSESS), and
understand how forecast skill changes with lead time.

Concept

In the previous section, synthetic observations were generated from the true climate state using a perfect-model experiment. Although these observations contain random errors, they provide our best estimate of the system at the observation time.

Forecasting begins by using an observation as the initial condition. Once initialized, the numerical model predicts the future evolution of the climate system without additional observational information.

Because the Lorenz model is chaotic, even a small error in the initial observation gradually amplifies with time. As a result, the forecast trajectory diverges from the true climate state, and forecast skill decreases as the forecast lead time increases.

Within a perfect-model framework, both the forecast and the true climate state are available, allowing forecast errors to be evaluated objectively. This provides an ideal environment for understanding the fundamental limits of predictability before introducing ensemble prediction and data assimilation.

Model / Equation

Model

The numerical model used in this chapter is identical to that used in the previous chapters. The Lorenz climate model, model parameters, and numerical integration remain unchanged. External forcing is set to zero,

AMP_X="0.0"
AMP_Y="0.0"
AMP_Z="0.0"

The only difference from the previous chapter is the experiment design. Instead of generating synthetic observations, we now use the observations to initialize forecasts and evaluate their predictive skill.

Forecast and Hindcast Experiments

A forecast experiment begins from an observed initial condition and predicts the future evolution of the climate system.

The forecast is initialized from the observation,

\[ \mathbf{X}_f(t_{\rm initial},0) = \mathbf{X}_{obs}(t_{\rm initial}), \]

where \(t_{\rm initial}\) denotes the forecast initialization time and the forecast lead time is initially zero.

As the model integrates forward,

\[ \mathbf{X}_f(t_{\rm initial},t_{\rm lead}), \]

describes the forecast at lead time \(t_{\rm lead}\).

In an operational forecast, the future climate state has not yet occurred, so forecast accuracy cannot be evaluated until observations become available.

In contrast, a hindcast experiment uses exactly the same forecast procedure, but the initialization time is selected from the past. Since the future evolution following that initialization has already occurred, the forecast can be compared directly with the known climate state.

Therefore, the only difference between a forecast and a hindcast is whether the future has already occurred.

The true climate state generated by the control simulation is

\[ \mathbf{X}_{true}(t), \]

where \(t\) is the model time.

For hindcast verification, it is convenient to reorganize the same trajectory using the initialization time and forecast lead time,

\[ t = t_{\rm initial} + t_{\rm lead}, \]

so that

\[ \boxed{ \mathbf{X}_{true}(t_{\rm initial},t_{\rm lead}) = \mathbf{X}_{true}(t_{\rm initial}+t_{\rm lead}) } \]

This equation does not define a new model. It simply expresses the same control simulation using forecast coordinates.

For example, a single control simulation,

Time
0--------------------------------------------------------------100

can be viewed as a collection of hindcasts,

Forecast 1

Initial time = 20

Lead time
0------1------2------3------4------...
│
▼
20-----21-----22-----23-----24-----...


Forecast 2

Initial time = 40

Lead time
0------1------2------3------4------...
│
▼
40-----41-----42-----43-----44-----...


Forecast 3

Initial time = 60

Lead time
0------1------2------3------4------...
│
▼
60-----61-----62-----63-----64-----...

Each hindcast begins from a different initialization time but is evaluated using the same forecast lead time. This common coordinate system allows forecast errors from many hindcasts to be combined and averaged to estimate predictive skill.

Predictive Skill

The forecast error is defined as

\[ \mathbf{E}(t_{\rm initial},t_{\rm lead}) = \mathbf{X}_f(t_{\rm initial},t_{\rm lead}) - \mathbf{X}_{true}(t_{\rm initial},t_{\rm lead}). \]

Forecast skill is evaluated by comparing many hindcasts initialized at different times. Throughout this section, the forecast lead time is fixed, while the initialization time varies.

Suppose that forecasts are initialized at

\[ t_{\rm initial}^{(1)}, t_{\rm initial}^{(2)}, \ldots, t_{\rm initial}^{(N)}, \]

where \(N\) is the total number of hindcast experiments.

The Root Mean Square Error (RMSE) at a given forecast lead time is

\[ {\rm RMSE}(t_{\rm lead}) = \sqrt{ \frac{1}{N} \sum_{i=1}^{N} \left[ X_f(t_{\rm initial}^{(i)},t_{\rm lead}) - X_{true}(t_{\rm initial}^{(i)},t_{\rm lead}) \right]^2 }. \]

Thus, RMSE is calculated by averaging the forecast errors over all hindcasts having the same forecast lead time.

The Anomaly Correlation Coefficient (ACC) at a given forecast lead time is

\[ {\rm ACC}(t_{\rm lead}) = \frac{ \displaystyle \frac{1}{N} \sum_{i=1}^{N} \left( X_f^{(i)} - \overline{X_f} \right) \left( X_{true}^{(i)} - \overline{X_{true}} \right) }{ \sigma_f \sigma_{true} }, \]

where

\[ X_f^{(i)} = X_f(t_{\rm initial}^{(i)},t_{\rm lead}), \]

\[ X_{true}^{(i)} = X_{true}(t_{\rm initial}^{(i)},t_{\rm lead}), \]

and

\[ \sigma_f \]

and

\[ \sigma_{true} \]

are the standard deviations of the forecast and true climate state computed across all initialization times.

ACC measures how well the forecasts reproduce the temporal variability of the true climate state across the collection of hindcasts.

Finally, the RMSE Skill Score is defined as

\[ {\rm RMSE\ Skill\ Score}(t_{\rm lead}) = 1- \frac{ {\rm RMSE}(t_{\rm lead}) } { {\rm RMSE}_{ref} }. \]

In this chapter,

\[ {\rm RMSE}_{ref} = \sigma_{true}, \]

where \(\sigma_{true}\) is the standard deviation of the control simulation. This reference represents the natural variability of the climate system and provides a simple benchmark for evaluating forecast skill.

A perfect forecast has

\[ {\rm RMSE}=0, \]

\[ {\rm ACC}=1, \]

and

\[ {\rm RMSE\ Skill\ Score}=1. \]

As forecast lead time increases, RMSE generally increases, while ACC and the RMSE Skill Score decrease because forecast errors grow through nonlinear dynamics.

Sample Code

The sample programs in this section reuse the Model_src-v3.0 framework developed in the previous exercises. The control simulation provides the true climate state, while multiple hindcast experiments are initialized from synthetic observations with prescribed observation error.

Download the source code archive,

Model_src-v3.0.tar.gz

and extract it using

tar -xzvf Model_src-v3.0.tar.gz

This creates the directory

Model_src-v3.0/
├── figures/
│   ├── fig-timeseries_spaghetti_simple.py
│   ├── fig-timeseries_hindcast_acc.py
│   └── fig-hindcast_skill-template.py
├── utils/
│   ├── get_observation.py
│   └── create_initial_obs.py
├── Lmodel_climate-v3.py
└── run_model_climate-v3.sh

The individual source files used in this section are available below.

File	Description
`run_model_climate-v3.sh`	Defines the experiment settings and executes the Lorenz climate model.
`Lmodel_climate-v3.py`	Performs numerical integration of the Lorenz climate model.
`get_observation.py`	Generates synthetic observations from the control simulation by adding random observation error.
`create_initial_obs.py`	Creates hindcast initial conditions from the synthetic observations at a prescribed initialization time.
`fig-timeseries_spaghetti_simple.py`	Plots the control simulation and multiple hindcast trajectories in a spaghetti diagram.
`fig-timeseries_hindcast_acc.py`	Visualizes the true and forecast states at selected lead times and displays the corresponding ACC values.
`fig-hindcast_skill-template.py`	Organizes the hindcast data by initialization time and lead time and calculates ACC, RMSE, and RMSE Skill Score for `X`, `Y`, and `Z`.

The hindcast experiments are initialized at

T_INITIAL = 55, 60, 65, 70, 75, 80, 85, and 90

and stored in separate directories:

HINDCAST/
├── 0055/
├── 0060/
├── 0065/
├── 0070/
├── 0075/
├── 0080/
├── 0085/
├── 0090/
└── run_hindcast.sh

The forecast data are reorganized into two-dimensional arrays such as

Xf(t_initial, t_lead)
Xtrue(t_initial, t_lead)

and similarly for Y and Z.

The first dimension represents the initialization experiment, while the second dimension represents forecast lead time. This data structure makes it straightforward to calculate forecast skill at each lead time across all hindcast experiments.

The following three metrics are evaluated:

ACC measures the linear correspondence between the forecast and true states.
RMSE measures the magnitude of the forecast error.
RMSE Skill Score compares the hindcast RMSE with a climatological reference forecast.

The RMSE Skill Score is interpreted as

RMSESS = 1    perfect forecast
RMSESS > 0    better than climatology
RMSESS = 0    equal to climatology
RMSESS < 0    worse than climatology

The final hindcast-skill figure contains nine panels showing ACC, RMSE, and RMSE Skill Score for X, Y, and Z as functions of forecast lead time.

Exercise

In this exercise, you will develop a simple forecast system using the perfect-model framework introduced in the previous chapter. The control simulation is regarded as the true climate state, while synthetic observations are used to initialize forecasts.

You will first generate the control simulation and synthetic observations. Next, you will perform forecast experiments initialized from the observations and compare the forecasts with the corresponding true climate states. Finally, you will extend the experiment to a hindcast framework and evaluate predictive skill using RMSE, ACC, and the RMSE Skill Score.

By the end of this exercise, you will understand how forecast skill is quantified and why forecast accuracy generally decreases as forecast lead time increases.

Step 1. Generate the Perfect-Model Dataset

In this chapter, we use the perfect-model framework developed in the previous section. First, perform a control simulation to generate the true climate state. Then, create synthetic observations from the control simulation.

Create the experiment directory,

mkdir Chapter8_3
cd Chapter8_3

Create directories for the control simulation and visualization,

mkdir CTL
mkdir Figure

The directory structure should look like

Chapter8_3/
├── CTL/
└── Figure/

Move to the control-simulation directory,

cd CTL

Copy the model source code,

cp ../../Model_src-v3.0/*.* .

Open run_model_climate-v3.sh and use

REST="False"

TSTR="0"
TEND="100"

AMP_X="0.0"
AMP_Y="0.0"
AMP_Z="0.0"

DEBUG_OUTPUT="False"

Run the control simulation,

./run_model_climate-v3.sh

or

bash run_model_climate-v3.sh

This simulation produces the true climate trajectory,

\[ \mathbf{X}_{true}(t), \]

which will be used for both forecast initialization and forecast verification.

Next, copy the observation-generation program,

cp ../../Model_src-v3.0/utils/get_observation.py .

Open get_observation.py and use

TSTR = 50
TEND = 100

OBS_DT = 5.0

OBS_TSTR = 55.0
OBS_TEND = 95.0

OERR = 2.0

outfile_nc = "observation.nc"

where

TSTR and TEND define the time range of the control simulation used to generate the synthetic observations.
OBS_DT specifies the observation interval.
OBS_TSTR and OBS_TEND define the observation period.
OERR specifies the standard deviation of the Gaussian observation error.

Important: The observation period must lie within the control simulation period. Therefore, the following condition must always be satisfied:

TSTR ≤ OBS_TSTR < OBS_TEND ≤ TEND

In this example,

Control simulation : T = 50–100
Observations       : T = 55–95

The observations are generated only from the available control simulation. If OBS_TSTR is earlier than TSTR or OBS_TEND is later than TEND, the program cannot generate observations outside the available model output.

Run the program,

python get_observation.py

This creates

observation.nc

which contains

observation times,
sampled true climate states, and
synthetic observations with random observation errors.

Next, move to the figure directory,

cd ../Figure

Copy the plotting program,

cp ../../Model_src-v3.0/figures/fig-timeseries_obs.py .

Open fig-timeseries_obs.py and confirm that

FDIR = "../CTL/out/"
OBSFILE = "../CTL/observation.nc"

Run the program,

python fig-timeseries_obs.py

The resulting figure should display

the continuous model trajectory,
the sampled true states, and
the synthetic observations.

These observations will serve as the initial conditions for the forecast experiments in the following exercises.

Questions

Why is the control simulation regarded as the true climate state in a perfect-model experiment?
Why are observations generated only every five time units instead of every model time step?
Why do the observations differ from the sampled true state?
Which observation times will be used to initialize the forecast experiments?

Step 2. Perform the First Hindcast Experiment

Next, perform the first hindcast experiment using the synthetic observation at model time 55 as the initial condition. The model will be integrated from

\[ t_{\rm initial}=55 \]

to

\[ t=60. \]

Thus, the maximum forecast lead time is

\[ t_{\rm lead}=5. \]

Return to the main experiment directory,

cd ..

Create the hindcast directory,

mkdir HINDCAST
cd HINDCAST

Create a directory for the hindcast initialized at model time 55,

mkdir 0055
cd 0055

The directory structure should now look like

Chapter8_3/
├── CTL/
├── Figure/
└── HINDCAST/
    └── 0055/

Copy the model source code,

cp ../../../Model_src-v3.0/*.* .

Copy the utility program used to create an initial condition from the synthetic observations,

cp ../../../Model_src-v3.0/utils/create_initial_obs.py .

The program reads

../../CTL/observation.nc

and extracts the synthetic observation at the prescribed initialization time.

Before running the utility program, modify create_initial_obs.py directly to set the default initialization time.

Open

vi create_initial_obs.py

and change the default value of the initialization time to

DEFAULT_T_INITIAL = 55

Save the file.

The program is now configured to create the restart file for model time 55 by default. Create the restart file,

python create_initial_obs.py

This creates

rest/restart.0055.nc

Verify the restart file,

ncdump rest/restart.0055.nc

The restart file contains the observed initial state,

\[ \mathbf{X}_{f}(55,0) = \mathbf{X}_{obs}(55). \]

Open run_model_climate-v3.sh and set

REST="True"

TSTR="55"
TEND="60"

AMP_X="0.0"
AMP_Y="0.0"
AMP_Z="0.0"

DEBUG_OUTPUT="False"

Confirm that the restart-file settings point to

rest/restart.0055.nc

Run the hindcast simulation,

./run_model_climate-v3.sh

or

bash run_model_climate-v3.sh

The hindcast begins from the noisy observation at model time 55 and predicts the model state through model time 60.

The relationship between model time and forecast lead time is

\[ t_{\rm lead} = t-t_{\rm initial}. \]

For this experiment,

Model time	Forecast lead time
55	0
56	1
57	2
58	3
59	4
60	5

The hindcast output represents

\[ \mathbf{X}_{f}(55,t_{\rm lead}), \]

while the corresponding true state is obtained from the control simulation,

\[ \mathbf{X}_{true}(55,t_{\rm lead}) = \mathbf{X}_{true}(55+t_{\rm lead}). \]

Visualize the First Hindcast

Return to the figure directory,

cd ../../Figure

Copy the control–perturbed comparison template and rename it for the hindcast initialized at model time 55,

cp ../../Model_src-v3.0/figures/fig-timeseries_CTL_PRD-template.py \
   fig-timeseries_CTL_PRD_0055.py

Open fig-timeseries_CTL_PRD_0055.py and modify the input directories and plotting periods,

# Model output directories
FDIR_CTL = "../CTL/out/"
FDIR_PRD = "../HINDCAST/0055/out/"

# Use a longer period for the control simulation
TSTR_CTL = 50
TEND_CTL = 65

# Use the hindcast period
TSTR_PRD = 55
TEND_PRD = 60

Update the figure labels, where necessary, so that the second simulation is identified as the hindcast rather than the perturbed simulation. For example, change

label="PRD"

to

label="Hindcast"

Run the plotting program,

python fig-timeseries_CTL_PRD_0055.py

The figure compares the continuous true climate trajectory from the control simulation with the hindcast initialized from the synthetic observation at model time 55.

At the initialization time, the hindcast and true trajectory are close but not identical because the synthetic observation contains random error. As the forecast lead time increases, the hindcast may gradually diverge from the true climate trajectory.

Questions

What initial condition is used for this hindcast?
Why does the hindcast not begin from the exact true state?
What is the difference between model time and forecast lead time?
What is the maximum forecast lead time in this experiment?
Why can this experiment be evaluated immediately as a hindcast?
How closely does the hindcast follow the control simulation during the five-unit forecast period?
Does the hindcast error appear to increase with forecast lead time?

Step 3: Automating the Hindcast Experiment with a Bash Script

In Step 2, we manually created a hindcast experiment initialized at T = 55. In this step, we automate the same workflow. We first repeat the experiment manually for T_INITIAL = 60 to understand each command. We then combine the commands into a reusable Bash script and verify the script by creating a new hindcast initialized at T_INITIAL = 65.

The directory structure at the beginning of this exercise should be:

Chapter8_3/
├── CTL/
├── Figure/
└── HINDCAST/
    └── 0055/

Move to the HINDCAST directory:

cd HINDCAST
pwd
ls

The output should include:

Step 3.1. Manually Create the `T_INITIAL = 60` Hindcast

In this section, we execute every command manually to understand the complete hindcast workflow.

Step 3.1.1. Create the `0060` Directory

mkdir 0060
ls

Expected output:

0055
0060

Step 3.1.2. Move into the Experiment Directory

cd 0060

Step 3.1.3. Copy the Model Source Files

cp ../../../Model_src-v3.0/*.* .

Confirm that run_model_climate-v3.sh has been copied.

ls

Step 3.1.4. Create the Initial Condition

Copy the initialization script:

cp ../../../Model_src-v3.0/utils/create_initial_obs.py .

Generate the restart file:

python create_initial_obs.py --t_initial=60

This program reads the control simulation and creates the initial condition for the hindcast initialized at T = 60.

Step 3.1.5. Create a Restart Run Script

Rather than editing the original run_model_climate-v3.sh, we will use the Unix sed command to generate a modified copy containing the restart settings.

The restart experiment should use

REST="True"
TSTR="60"
TEND="65"

Generate the new run script:

sed \
    -e 's/^REST=.*/REST="True"/' \
    -e 's/^TSTR=.*/TSTR="60"/' \
    -e 's/^TEND=.*/TEND="65"/' \
    run_model_climate-v3.sh > run_model_climate-v3_0060.sh

The original file remains unchanged, while a new script (run_model_climate-v3_0060.sh) is created specifically for this hindcast experiment.

Step 3.1.6. Verify the Model Settings

grep -E '^(REST|TSTR|TEND)=' run_model_climate-v3_0060.sh

Expected output:

REST="True"
TSTR="60"
TEND="65"

Step 3.1.7. Run the Hindcast

bash run_model_climate-v3_0060.sh

After the model finishes,

ls out

The out directory should contain the model output for T = 60--65.

Step 3.2. Create an Automated Bash Script

After confirming that the manual experiment works correctly, we can automate the entire workflow using a Bash script. The script will create a new experiment directory, generate the initial condition, configure the model run, and execute the hindcast automatically.

Return to the HINDCAST directory:

cd ..
pwd

Create the script:

vi run_hindcast.sh

Paste the following code.

#!/bin/bash

# ============================================================
# Automated Hindcast Experiment
#
# Usage:
#   ./run_hindcast.sh T_INITIAL
#
# Example:
#   ./run_hindcast.sh 65
# ============================================================

set -e

if [ $# -ne 1 ]; then
    echo "Usage: $0 T_INITIAL"
    exit 1
fi

T_INITIAL=$1
T_END=$((T_INITIAL + 5))

EXP_NAME=$(printf "%04d" "${T_INITIAL}")
RUN_SCRIPT="run_model_climate-v3_${EXP_NAME}.sh"

echo "========================================"
echo "Create Hindcast Experiment"
echo "========================================"
echo "T_INITIAL : ${T_INITIAL}"
echo "T_END     : ${T_END}"
echo "Directory : ${EXP_NAME}"
echo "========================================"

mkdir -p "${EXP_NAME}"
cd "${EXP_NAME}"

cp ../../../Model_src-v3.0/*.* .
cp ../../../Model_src-v3.0/utils/create_initial_obs.py .

python create_initial_obs.py --t_initial="${T_INITIAL}"

sed \
    -e 's/^REST=.*/REST="True"/' \
    -e "s/^TSTR=.*/TSTR=\"${T_INITIAL}\"/" \
    -e "s/^TEND=.*/TEND=\"${T_END}\"/" \
    run_model_climate-v3.sh > "${RUN_SCRIPT}"

echo "----------------------------------------"
echo "Model Settings"
echo "----------------------------------------"
grep -E '^(REST|TSTR|TEND)=' "${RUN_SCRIPT}"
echo "----------------------------------------"

bash "${RUN_SCRIPT}"

cd ..

echo "========================================"
echo "Hindcast experiment ${EXP_NAME} completed."
echo "========================================"

Save the file and make it executable:

chmod +x run_hindcast.sh

Step 3.3. Test the Bash Script with `T_INITIAL = 65`

At this point, the 0055 and 0060 experiments already exist. We will verify that the Bash script works by creating a new experiment initialized at T_INITIAL = 65.

Run

./run_hindcast.sh 65

The script will automatically:

Create the experiment directory 0065.
Copy the model source files.
Generate the initial condition at T = 65.
Create run_model_climate-v3_0065.sh.
Set REST=True, TSTR=65, and TEND=70.
Execute the hindcast simulation.

Confirm the new directory:

ls

Expected output:

0055
0060
0065
run_hindcast.sh

Verify the generated run script:

grep -E '^(REST|TSTR|TEND)=' 0065/run_model_climate-v3_0065.sh

Expected output:

REST="True"
TSTR="65"
TEND="70"

Summary

During this exercise, the Bash script automates the following workflow:

CTL simulation
      │
      ▼
create_initial_obs.py
      │
      ▼
Restart file
      │
      ▼
sed
      │
      ▼
run_model_climate-v3_0065.sh
      │
      ▼
bash
      │
      ▼
Hindcast simulation

Step 4: Visualizing Multiple Hindcast Experiments

In Step 3, we created three hindcast experiments initialized at:

T_INITIAL = 55
T_INITIAL = 60
T_INITIAL = 65

The corresponding experiment directories are:

HINDCAST/
├── 0055/
├── 0060/
└── 0065/

In Step 4, we will first visualize these three hindcast experiments using a spaghetti plot. We will then create additional hindcasts initialized at T = 70, 75, 80, 85, and 90 and update the spaghetti plot to include all experiments.

The correct spelling is spaghetti, so the figure script is named:

fig-timeseries_spaghetti_simple.py

The template script is archived in:

Model_src-v3.0/figures/

Step 4.1. Copy the Spaghetti-Plot Script

Move to the Figure directory:

cd Figure

Copy the archived template script:

cp ../../../Model_src-v3.0/figures/fig-timeseries_spaghetti_simple.py .

Confirm that the script was copied:

ls

The output should include:

fig-timeseries_spaghetti_simple.py

Step 4.2. Create the Initial Spaghetti Plot

The archived script is configured to read:

../CTL/out/
../HINDCAST/0055/out/
../HINDCAST/0060/out/
../HINDCAST/0065/out/

Run the script:

python fig-timeseries_spaghetti_simple.py

The output figure is:

fig-timeseries_spaghetti_simple.png

In this figure:

the black curve shows the control simulation;
the orange curves show the hindcast experiments; and
the circles indicate the beginning of each hindcast.

Open the figure and confirm that the three hindcast experiments are displayed.

Step 4.3. Create Additional Hindcast Experiments

Return to the HINDCAST directory:

cd ../HINDCAST

Use the automated Bash script developed in Step 3 to create additional hindcast experiments.

Create a hindcast initialized at T = 70:

./run_hindcast.sh 70

Create a hindcast initialized at T = 75:

./run_hindcast.sh 75

Create a hindcast initialized at T = 80:

./run_hindcast.sh 80

Create a hindcast initialized at T = 85:

./run_hindcast.sh 85

Create a hindcast initialized at T = 90:

./run_hindcast.sh 90

Each hindcast is integrated for five model-time units. Therefore, the experiment periods are:

0070: T = 70–75
0075: T = 75–80
0080: T = 80–85
0085: T = 85–90
0090: T = 90–95

Confirm the experiment directories:

ls

The output should include:

0055
0060
0065
0070
0075
0080
0085
0090
run_hindcast.sh

Step 4.4. Update the Spaghetti-Plot Script

Return to the Figure directory:

cd ../Figure

Open the script:

vi fig-timeseries_spaghetti_simple.py

Change the ending time displayed in the figure from:

TEND_CTL = 75

to:

TEND_CTL = 95

Next, update the list of hindcast experiments.

Replace:

HINDCAST_EXPERIMENTS = [
    "0055",
    "0060",
    "0065",
]

with:

HINDCAST_EXPERIMENTS = [
    "0055",
    "0060",
    "0065",
    "0070",
    "0075",
    "0080",
    "0085",
    "0090",
]

Save the revised script and run it again:

python fig-timeseries_spaghetti_simple.py

The updated figure will contain eight hindcast trajectories initialized every five model-time units from T = 55 to T = 90.

The output file remains:

fig-timeseries_spaghetti_simple.png

Step 4.5. Examine the Spaghetti Diagram

Compare each hindcast trajectory with the control simulation.

Consider the following questions:

Does each hindcast begin near the control trajectory?
How quickly does each hindcast separate from the control simulation?
Do the forecast errors grow similarly for all initialization times?
Are some initial conditions more predictable than others?
How does the nonlinear evolution of the Lorenz model affect forecast predictability?

The spaghetti diagram shows that forecast error growth depends not only on forecast lead time but also on the model state at initialization. Some hindcasts may remain close to the control trajectory, whereas others may diverge rapidly.

Step 5: Evaluating Hindcast Skill

In Step 4, we created multiple hindcast experiments initialized every five model-time units:

In Step 5, we evaluate the predictive skill of these hindcasts using three metrics:

Anomaly Correlation Coefficient (ACC)
Root-Mean-Square Error (RMSE)
RMSE Skill Score (RMSESS)

Two template scripts are provided:

fig-timeseries_hindcast_acc.py
fig-hindcast_skill-template.py

The first script illustrates how ACC is calculated at selected lead times. The second script organizes all hindcast data by initialization time and lead time and then calculates ACC, RMSE, and RMSESS for X, Y, and Z.

The template scripts are archived in:

Model_src-v3.0/figures/

Step 5.1. Copy the Hindcast-Skill Scripts

Move to the Figure directory:

cd Figure

Copy the ACC visualization script:

cp ../../../Model_src-v3.0/figures/fig-timeseries_hindcast_acc.py .

Copy the full hindcast-skill template:

cp ../../../Model_src-v3.0/figures/fig-hindcast_skill-template.py .

Confirm that both scripts were copied:

ls

The output should include:

fig-timeseries_hindcast_acc.py
fig-hindcast_skill-template.py

Step 5.2. Visualize ACC at Selected Lead Times

First, run:

python fig-timeseries_hindcast_acc.py

The output figure is:

fig-timeseries_hindcast_acc.png

This figure contains five panels corresponding to:

Lead time = 0
Lead time = 1
Lead time = 2
Lead time = 3
Lead time = 4

Each panel shows:

the control simulation as the true trajectory;
all hindcast trajectories;
the true state at the selected verification time;
the corresponding forecast state;
a line connecting the true and forecast states; and
the ACC value calculated across all hindcast experiments.

ACC is calculated using:

np.corrcoef(true_values, forecast_values)[0, 1]

The correlation is calculated across all initialization experiments at one common lead time.

For example, at lead time 2, the comparison includes:

Experiment 0055: truth and forecast at T = 57
Experiment 0060: truth and forecast at T = 62
Experiment 0065: truth and forecast at T = 67
...
Experiment 0090: truth and forecast at T = 92

Examine how the forecast trajectories separate from the control simulation as lead time increases.

Step 5.3. Organize Hindcast Data by Initialization Time and Lead Time

The second template reorganizes each model variable into two-dimensional arrays:

Xf[t_initial, t_lead]
Xtrue[t_initial, t_lead]

and similarly:

Yf[t_initial, t_lead]
Ytrue[t_initial, t_lead]

Zf[t_initial, t_lead]
Ztrue[t_initial, t_lead]

The first dimension represents the hindcast initialization experiment:

0055, 0060, 0065, ..., 0090

The second dimension represents forecast lead time.

For example:

Xf[0, :]

contains the complete X forecast trajectory initialized at T = 55.

Xf[:, 100]

contains the X forecasts from all initialization experiments at the same lead-time index.

This structure makes it straightforward to calculate skill across all hindcast experiments at each lead time.

Step 5.4. Define the Lead-Time Coordinate

The lead-time coordinate is derived from the control simulation so that it has the same temporal resolution as the model output.

The template uses:

FORECAST_LENGTH = 5.0

lead_time_ctl = time_ctl.sel(
    time=slice(
        TSTR_CTL,
        TSTR_CTL + FORECAST_LENGTH
    )
)

LEAD_TIMES = lead_time_ctl.values[:-1] - TSTR_CTL

The final endpoint is excluded because each hindcast contains output from the initialization time up to, but not including, the exact ending time.

For example, if the model output interval is 0.01, the lead times are approximately:

0.00, 0.01, 0.02, ..., 4.99

Check the lead-time information printed by the script:

Number of lead times
First lead time
Last lead time
Lead-time interval

Step 5.5. Calculate ACC

At each lead time, ACC is calculated across all hindcast initialization experiments.

For X, the calculation is:

ACC_X[ilead] = np.corrcoef(
    Xtrue[:, ilead],
    Xf[:, ilead]
)[0, 1]

The same calculation is applied to Y and Z.

ACC measures the similarity between forecast and true variations across initialization cases.

ACC = 1    perfect linear agreement
ACC = 0    no linear relationship
ACC = -1   perfect inverse relationship

Step 5.6. Calculate RMSE

RMSE measures the magnitude of forecast error.

For X, the calculation is:

RMSE_X[ilead] = np.sqrt(
    np.mean(
        (Xf[:, ilead] - Xtrue[:, ilead])**2
    )
)

The same calculation is applied to Y and Z.

A smaller RMSE indicates a more accurate forecast.

RMSE = 0    perfect forecast
larger RMSE greater forecast error

Step 5.7. Calculate the RMSE Skill Score

The RMSE Skill Score compares the hindcast RMSE against the RMSE of a climatological reference forecast.

The climatological reference is calculated from the control simulation:

Xclim = float(Xctl.mean())
Yclim = float(Yctl.mean())
Zclim = float(Zctl.mean())

For X, the climatological RMSE is:

RMSE_CLIM_X[ilead] = np.sqrt(
    np.mean(
        (Xclim - Xtrue[:, ilead])**2
    )
)

The RMSE Skill Score is:

RMSESS_X[ilead] = (
    1.0 - RMSE_X[ilead] / RMSE_CLIM_X[ilead]
)

The same calculation is applied to Y and Z.

Interpret the RMSE Skill Score as:

RMSESS = 1    perfect forecast
RMSESS > 0    better than climatology
RMSESS = 0    equal to climatology
RMSESS < 0    worse than climatology

Step 5.8. Run the Hindcast-Skill Template

Run:

python fig-hindcast_skill-template.py

The output figure is:

fig-hindcast_skill.png

The script removes _template from the output filename automatically.

The figure contains nine panels arranged as:

            ACC        RMSE       RMSE Skill Score

X            X           X                X

Y            Y           Y                Y

Z            Z           Z                Z

All panels use forecast lead time on the x-axis.

The skill values are calculated at the full temporal resolution of the model output, but only major lead-time ticks are displayed.

Step 5.9. Examine the Hindcast-Skill Figure

Examine how predictive skill changes with forecast lead time.

Consider the following questions:

Does ACC decrease as lead time increases?
Does RMSE increase as lead time increases?
At what lead time does RMSESS approach zero?
Are the predictability limits similar for X, Y, and Z?
Does one variable lose predictive skill more rapidly than the others?
Are there short lead-time fluctuations in the skill curves?
How does nonlinear error growth affect the three skill metrics differently?

ACC and RMSE describe different aspects of forecast quality. ACC measures the correspondence between forecast and true variations, whereas RMSE measures the magnitude of forecast error. RMSESS indicates whether the forecast improves upon a climatological reference.

Key Points

After completing this chapter, you should understand the following concepts.

A hindcast experiment evaluates model predictability by initializing forecasts from known model states.
Automating repetitive experiments with a Bash script greatly improves efficiency and reproducibility.
Forecast trajectories initialized from nearby states diverge because of the chaotic nature of the Lorenz system.
Organizing forecast data by initialization time and forecast lead time simplifies forecast verification.
ACC measures how well forecasts reproduce the temporal variations of the true state.
RMSE measures the magnitude of forecast errors.
RMSE Skill Score compares the forecast with a climatological reference and indicates whether the forecast provides additional predictive skill.
Forecast skill generally decreases as forecast lead time increases because initial-condition errors amplify in a chaotic system.
Observation error influences the quality of the initial conditions and therefore affects forecast skill.

Homework

In this homework, you will examine how reducing observation error influences hindcast skill.

The original exercise uses

OERR = 2.0

You will repeat the complete hindcast experiment using

OERR = 1.0

and compare the results.

Task

Open get_observation.py and change

OERR = 2.0

to

OERR = 1.0

Keep the same observation period and observation interval used in the original experiment.

Generate the new synthetic observation dataset:

python get_observation.py

Next, rerun the hindcast experiments initialized at

T_INITIAL = 55, 60, 65, 70, 75, 80, 85, and 90

Use the Bash script developed in Step 3:

./run_hindcast.sh 55
./run_hindcast.sh 60
./run_hindcast.sh 65
./run_hindcast.sh 70
./run_hindcast.sh 75
./run_hindcast.sh 80
./run_hindcast.sh 85
./run_hindcast.sh 90

After all hindcasts have been completed, move to the Figure directory and run

python fig-hindcast_skill-template.py

This creates

fig-hindcast_skill.png

containing ACC, RMSE, and RMSE Skill Score for X, Y, and Z.

Questions

Compare the new experiment with OERR = 1.0 against the original experiment with OERR = 2.0.

Answer the following questions briefly.

How does reducing OERR from 2.0 to 1.0 affect ACC?
How does reducing OERR affect RMSE?
How does reducing OERR affect RMSE Skill Score?
Does the OERR = 1.0 experiment remain skillful for a longer forecast lead time?
Are the effects of reducing observation error similar for X, Y, and Z?
Why does reducing observation error improve the quality of the hindcast initial conditions?
How does the improvement in the initial conditions influence forecast error growth in the Lorenz system?

Submission

Submit the following figure:

fig-hindcast_skill.png generated using OERR = 1.0

Also submit one short paragraph comparing the hindcast-skill results obtained with

OERR = 1.0

and

OERR = 2.0

Your paragraph should discuss changes in ACC, RMSE, and RMSE Skill Score and explain how the observation error influences the quality of the initial conditions and subsequent forecast skill.

Chapter 8.3 Forecast Experiment and Predictive Skill

Learning Objectives

Concept

Model / Equation

Model

Forecast and Hindcast Experiments

Predictive Skill

Sample Code

Exercise

Step 1. Generate the Perfect-Model Dataset

Step 2. Perform the First Hindcast Experiment

Visualize the First Hindcast

Step 3: Automating the Hindcast Experiment with a Bash Script

Step 3.1. Manually Create the T_INITIAL = 60 Hindcast

Step 3.1.1. Create the 0060 Directory

Step 3.1.2. Move into the Experiment Directory

Step 3.1.3. Copy the Model Source Files

Step 3.1.4. Create the Initial Condition

Step 3.1.5. Create a Restart Run Script

Step 3.1.6. Verify the Model Settings

Step 3.1.7. Run the Hindcast

Step 3.2. Create an Automated Bash Script

Step 3.3. Test the Bash Script with T_INITIAL = 65

Summary

Step 4: Visualizing Multiple Hindcast Experiments

Step 4.1. Copy the Spaghetti-Plot Script

Step 4.2. Create the Initial Spaghetti Plot

Step 4.3. Create Additional Hindcast Experiments

Step 4.4. Update the Spaghetti-Plot Script

Step 4.5. Examine the Spaghetti Diagram

Step 5: Evaluating Hindcast Skill

Step 5.1. Copy the Hindcast-Skill Scripts

Step 5.2. Visualize ACC at Selected Lead Times

Step 5.3. Organize Hindcast Data by Initialization Time and Lead Time

Step 5.4. Define the Lead-Time Coordinate

Step 5.5. Calculate ACC

Step 5.6. Calculate RMSE

Step 5.7. Calculate the RMSE Skill Score

Step 5.8. Run the Hindcast-Skill Template

Step 5.9. Examine the Hindcast-Skill Figure

Key Points

Homework

Task

Questions

Submission

Step 3.1. Manually Create the `T_INITIAL = 60` Hindcast

Step 3.1.1. Create the `0060` Directory

Step 3.3. Test the Bash Script with `T_INITIAL = 65`