In a time series regression model I am working on, our intervention happened just several months after an big influential external event. The full time series model, therefore, underestimates the external event, but puts more weights on our intervention effect. However, we know that this is not quite correct from other researchs and domain experts. So the current remedy is to give an estimated effect to the external event associated with a given distribution, say the Normal distribution. After simulating for $N$ times, we plug this simulated number as a fixed value in the model, the rest of the coefficients are still estimated from the R software (using "arima" function). Therefore they result in $N$ models, or less.
Based on this approach, the distribution of some coefficients are not normally distributed. My question is:
Should I still choose the model based on the model selection criteria, or should I summarize each coefficient as a distribution?
There is a parameter originally not statistically significant, but with the expected sign. I kept it there as a confounding because it could change the estimate of the intervention parameter by 10%. If we summarize coefficients as distributions, the standard deviation of this distribution is much smaller than the standard error of the coefficient in the full model. And the 95% confidence range does not contain zero as well. Therefore how should I treat it for the reporting purpose, which goes back to Question 1, how should I summarize the models due to this method?
It will be really appreciated if someone can provide some suggestion.