On checking, it does seem that u_err and d_err are simply normal distributions, so I can sample in the proposed way. I may have been misled earlier because I was using a small number of samples and u_err and d_err did not seem to be normal.
The other thing which surprises me is how much better the fit is to training data compared with using a simple lagged variable - using a structural time series looks too good to be true - but the forecasting does have reasonable uncertainty.
Correction - on checking again, d_err is certainly not normal. The mean over time for each sample is negative.
so my impression is that for each sample there is a normal distribution over time with zero mean, but each sample is sampled from a distribution and accepted/rejected according to the likelihood, or something like that. Why is
sd(apply(d_err, 1, mean))
If this low value is correct, then we need to account for it in the forecasting?
On further investigation, looking at the diagnostics file, it turn out that each time step is treated as a separate parameter, so if you have 2000 time steps there are 2000 independent u_err parameters and 2000 independent d_err parameters. I am no longer surprised I get such a good fit to historical data. This topic needs a lot more further investigation.