# How to Determine whether Simulation Draws are Correct

I have implemented algorithm 1 and 2 in this paper http://www.lse.ac.uk/statistics/documents/researchreport61.pdf for the analysis/simulation hidden states for some time series. The reason why I am attempting to do this is so that I can automate the pre-match betting on sporting events and this is part of a much larger code base.

Now, I have implemented the algorithms outlined in the paper and the output of these is of course stochastic. For algorithm two I have the following pretty picture for a 1D observation vector time-series and a 1D state vector

1. Nile = actual data
2. alpha_hat = smoothed state
3. alpha_tilde = simulated smoothed state

Ok, the problem I have is that qualitatively I can see that the simulated smoothed state "alpha_tilde" is behaving "reasonably" and roughly following "alpha_hat" and indeed the observations, which I want it to do. However, when this is re-run I get a different, but also (via my eye) reasonable output.

I have some experience in analysis of periodicity of time-series (FFT, Wavelets, structure functions et al.), but what I want here is to determine quantitatively if my simulated state is correct. This is difficult because of it's random nature and it is not as simple as merely applying univariate statistics to the two series. I thought of using a variogram but I don't see how I can use it to compare the simulated series with the expected series and get a meaningful correlation?

Any advice on how to compare two time-series for correlation would be appreciated.

I don't know how you could use cross-correlation to quantitatively check the correctness of the implementation. Below I give you some other ideas.

One possibility (probably the most straightforward) to test your implementation of the simulation smoother is trying to reproduce the results from some application that use it for the model you are interested in.

The simulation smoother can be used for Bayesian estimation of state space models using the Gibbs sampler. Another possibility to test your code is therefore to generate several series from a known model and check if the average parameter estimates based on the Gibbs sampler are close to the true values (e.g., you could check if the true values lie within the 95% confidence intervals in 95% of the simulated series). This would require the additional work of implementing the Gibbs sampler, but for the local level model it is not that hard (at least for the purposes of getting some benchmark results).

Edit

If you are looking for some software to get benchmark results, the function simulateSSM in the R package KFAS implements the simulation smoother. In order to replicate results you may set the seed for the random number generator (and use the same algorithm to generate the random values from the Gaussian distribution). Alternatively you may dig into the source code and modify the function to pass some values or to print the values that are generated by the random number generator so that you can use them in your function to compare results. The first approach seems the most straightforward, but you will have to make sure that a given seed generates the same draws in R and in the software you are using.

Be also aware that the simulation smoother depends on the output from the Kalman filter (the prediction error, its variance and the Kalman gain) so you will have to check that the same input is passed to the simulation smoother (differences for example in the initialization of the Kalman filter may lead to changes in the output from the simulation smoother).

• Thank for your response, it is most appreciated. This is part of the code that is very important and it is also the only part where testing it is very difficult [to quantify]. The problem with the first suggesting (which I have thought about in the past) is A. finding software that can perform such operations atomically (just the simulation) and B. The random nature of the output. I think I will go with your second suggestion for now: simulate N series and take the mean for each point; see whether these mean values lie within the 95% confidence interval. Jan 4, 2015 at 10:32
• PS. I am currently calculating the confidence intervals using the standard formulation (same one that is on Wikipedia). This gives the same values above and below the series. Should i be doing this differently for State Space Time Series? Jan 4, 2015 at 10:34
• @Killercam I have edited the question as response to your first comment. As regards the confidence intervals, I think you could simply take the 0.025 and 0.975 quantiles from all the estimates (one for each simulated series) obtained for the variance parameter in the level component. For each simulated series your should record whether the true value of the parameter lied within that interval. I would expect it will do so in 95% of the cases. Jan 4, 2015 at 11:14
• Nice, thanks for this. I have KFAS so will take a look. However, I am not overly familiar with R; couple this with the number of variables that could effect the outcome of simulateSSM AND the random nature of the simulation process I am going to try options 2... Thanks very much for your time, it is most appreciated. Jan 4, 2015 at 11:49
• Finally, I wondered if you could clear something up. I am seeing graphs all over the place, where the confidence interval for a one-dimensional data set seems to change for each point in the plot; what type of confidence interval is this and should I be using it? Currently I have implemented the standard version which outputs a +/- value for the entire data set/series. Jan 4, 2015 at 12:34