diff --git a/probabilistic-analysis/probabilistic_assessment.py b/probabilistic-analysis/probabilistic_assessment.py
index 98a97d1..903f164 100644
--- a/probabilistic-analysis/probabilistic_assessment.py
+++ b/probabilistic-analysis/probabilistic_assessment.py
@@ -475,17 +475,18 @@ def process(beach_name, beach_scenario, n_runs, start_year, end_year,
                                               xp=profile_volume[valid_idx],
                                               fp=profile_chainage[valid_idx])
 
-            fig, ax = plt.subplots(9,
-                                   len(output_years),
-                                   figsize=(16, 24),
-                                   sharey='row')
-
             # Check whether to save probabilistic diagnostics
+            output_diagnostics = False
             for _, bp in pd.DataFrame(diagnostics).iterrows():
                 if ((str(prof['block']) == str(bp['block']))
                         and (prof['profile'] == bp['profile'])):
                     output_diagnostics = True
 
+                    fig, ax = plt.subplots(9,
+                                           len(output_years),
+                                           figsize=(16, 24),
+                                           sharey='row')
+
             # Loop through years
             pbar_year = tqdm(output_years, leave=False)
             for j, year in enumerate(pbar_year):
@@ -682,6 +683,7 @@ def process(beach_name, beach_scenario, n_runs, start_year, end_year,
                         'diagnostics',
                         '{} {} {}.csv'.format(beach_scenario, year,
                                               profile_type))
+                    dump_df = dump_df[::100]  # Only output every 100th row
                     dump_df.to_csv(csv_name, float_format='%g')
 
                     for i, c in enumerate(dump_df.columns[3:]):
@@ -714,6 +716,115 @@ def process(beach_name, beach_scenario, n_runs, start_year, end_year,
                 plt.savefig(figname, bbox_inches='tight', dpi=300)
                 plt.close(fig)
 
+                # Plot time series figure
+                fig, ax = plt.subplots(4,
+                                       2,
+                                       figsize=(12, 16),
+                                       sharey='row',
+                                       gridspec_kw={
+                                           'wspace': 0.05,
+                                           'width_ratios': [3, 1]
+                                       })
+
+                ax[0, 0].plot(years, slr[:, :100], c='#aaaaaa', lw=0.2)
+                ax[0, 0].plot(years,
+                              slr[:, 1],
+                              c='C0',
+                              label='Sample simulation')
+                ax[0, 1].hist(slr[-1, :],
+                              bins=100,
+                              fc='#cccccc',
+                              ec='#aaaaaa',
+                              orientation='horizontal')
+
+                ax[1, 0].plot(years, (slr * bf)[:, :100], c='#aaaaaa', lw=0.2)
+                ax[1, 0].plot(years, (slr * bf)[:, 1], c='C0')
+                ax[1, 1].hist((slr * bf)[-1, :],
+                              bins=100,
+                              fc='#cccccc',
+                              ec='#aaaaaa',
+                              orientation='horizontal')
+
+                ax[2, 0].plot(years, ur[:, :100], c='#aaaaaa', lw=0.2)
+                ax[2, 0].plot(years, ur[:, 1], c='C0')
+                ax[2, 1].hist(ur[-1, :],
+                              bins=100,
+                              fc='#cccccc',
+                              ec='#aaaaaa',
+                              orientation='horizontal')
+
+                ax[3, 0].plot(years, r[:, :100], c='#aaaaaa', lw=0.2)
+                ax[3, 0].plot(years, r[:, 1], c='C0', zorder=3)
+                ax[3, 1].hist(r[-1, :],
+                              bins=100,
+                              fc='#cccccc',
+                              ec='#aaaaaa',
+                              orientation='horizontal')
+
+                baseline = r[:, 1]
+                sdd = storm_demand_dist[:, 1]
+                for i in range(len(slr)):
+                    pe = [
+                        matplotlib.patheffects.Stroke(linewidth=5,
+                                                      foreground='#ffffff',
+                                                      capstyle='butt'),
+                        matplotlib.patheffects.Normal()
+                    ]
+
+                    ax[3, 0].plot([years[i], years[i]],
+                                  [baseline[i], baseline[i] + sdd[i]],
+                                  c='C0',
+                                  path_effects=pe)
+
+                # Maximum recession encountered
+                r_max = (baseline + sdd).max()
+                ax[3, 0].axhline(y=r_max, c='C3', linestyle=':')
+
+                i = len(years) - 1
+                for a in ax[:, 0]:
+                    a.set_xlim(right=years[-1])
+
+                # Add line at zero
+                for a in ax[:-1, 0]:
+                    a.spines['bottom'].set_position(('data', 0))
+                    a.set_xticklabels([])
+
+                for a in ax[:, 1]:
+                    a.xaxis.set_visible(False)
+                    a.spines['bottom'].set_visible(False)
+
+                ax[3, 0].annotate(
+                    'Most eroded beach state encountered in planning period',
+                    (years[0], r_max),
+                    xytext=(0, 20),
+                    textcoords='offset pixels')
+
+                ax[0, 0].legend()
+
+                ax[0, 0].set_title((f'Probabilistic trajectories\n'
+                                    f'(first 100 out of {n_runs:,} runs)'))
+                ax[0, 1].set_title(
+                    f'Probability\ndistribution\nin year {years[i]}')
+                ax[0, 0].set_ylabel('Sea level (m)', labelpad=10)
+                ax[1, 0].set_ylabel('Bruun recession (m)', labelpad=10)
+                ax[2, 0].set_ylabel('Underlying recession (m)', labelpad=10)
+                ax[3, 0].set_ylabel('Shoreline displacement (m)', labelpad=10)
+
+                ax[3, 0].set_title(('Bruun recession'
+                                    '+ underlying recession'
+                                    '+ storm demand'),
+                                   y=0.9)
+
+                for a in ax.ravel():
+                    a.spines['top'].set_visible(False)
+                    a.spines['right'].set_visible(False)
+
+                figname = os.path.join(
+                    'diagnostics',
+                    f'{beach_scenario} {profile_type} timeseries.png')
+                plt.savefig(figname, bbox_inches='tight', dpi=300)
+                plt.close(fig)
+
 
 def main():