Bootstrapping residuals: Am I doing it right? First of all: 
From what I understood, bootstrapping residuals works as follows:


*

*Fit model to data 

*Calculate the residuals

*Resample the residuals and add them to 1.

*Fit model to new dataset from 3.

*Repeat n times, but always add the resampled residuals to the fit
from 1.


Is that correct so far?

What I want to do is something slightly different:
I want to estimate parameter and prediction uncertainty for an algorithm that estimates some environmental variable.
What I have is a error-free time-series (from a simulation) of that variable, x_true, to which I add some noise, x_noise, in order to generate a synthetic dataset x.
I then try to find optimal parameters by fitting my algorithm with the sum of squares sum((x_estimate - x_true)^2) (! not x_estimate - x !) as an objective function. In order to see how my algorithm performs and to create samples of my parameters' distributions, I want to resample x_noise, add it to x_true, fit my model again, rinse and repeat. Is that a valid approach to assess parameter uncertainty? Can I interpret the fits to the bootstrapped datasets as prediction uncertainty, or do I have to follow the procedure I posted above?
/edit: I think I haven't really made clear what my model does. Think of it as essentially something like a de-noising method. It's not a predictive model, it's an algorithm that tries to extract the underlying signal of a noisy time-series of environmental data.
/edit^2: For the MATLAB-Users out there, I wrote down some quick & dirty linear regression example of what I mean.
This is what I believe "ordinary" bootstrapping of residuals is (please correct me if I'm wrong): http://pastebin.com/C0CJp3d1
This is what I want to do: http://pastebin.com/mbapsz4c
 A: To see how an algorithm performs in terms of predictive accuracy/mean squared error, you probably need the Efron-Gong "optimism" bootstrap.  This is implemented for easy use in the R rms package.  See its functions ols, validate.ols, calibrate.
A: Here is the general (semi-parametric-bootstrap) algorithm in more detail:
$\text{B}$ = number of bootstraps
the model:
$y = x\beta + \epsilon$
let $\hat{\epsilon}$ be the residuals


*

*Run the regression and obtain the estimator(s) $\hat\beta$ and residuals $\hat\epsilon$.

*Resample the residuals with replacement and obtain the bootstrapped residual vector $\hat\epsilon_\text{B}$.

*Obtain the bootstrapped dependent variable by multiplying the estimator(s) from (1) with the original regressors and adding the bootstrapped residual: $y_\text{B} = x\hat\beta + \hat\epsilon_\text{B}$.

*Run the regression with the bootstrapped dependent variables and the original regressors, this gives the bootstrapped estimator, i.e. regress $y_B$ on $x$, this gives $\hat\beta_\text{B}$.

*Repeat the procedure $\text{B}$-times by going back to (2).

A: I'm not sure that my understanding is correct. But here is my suggestion to modify your code ("ordinary bootstrapping of residuals", lines 28-34) into:
for i = 2:n_boot  
x_res_boot = x_residuals( randi(n_data,n_data,1) );  
x_boot = x_res_boot+ x_best_fit;  
p_est(:, i) = polyfit( t, x_boot, 1 );  
x_best_fit2 = polyval( p_est(:, i), t );  
x_residuals = x_best_fit2 - x_boot;
x_best_fit=x_best_fit2;
end  

The idea is that each time you are using residuals not from the first run, but from the previous bootstrap fit. As for me, all other seems to be valid.
This is revised version that has been checked in MATLAB. Two errors have been fixed.
A: Spatiotemporal application
I applied @Sweetbabyjesus's answer. Unfortunately my model is more complex than an OLS regression (it's a distributed-lag 'hhh4 framework' spatiotemporal model from surveillance R package) causing instability when making such data changes, resulting in non-convergence messages.
If this happens to you then a second-best solution is producing jackknifed estimates by fitting the model to your dataset with one data point missing and repeating for $N$ to get a distribution of model estimates; in a spatiotemporal model context I interpret this as subtracting 1 from the response variable $y$ e.g. disease case counts while avoiding negative counts, so it is a trimmed subtraction function. Since I cannot fit the entire model if I make a single space x time cell go missing i.e. NA, subtracting by 1 count was the only option.
The mean of jackknifed estimates was very close to the model point estimates, however the jackknifed sd estimates underestimated the model sd by 1-3 orders of magnitude! Generally speak jackknifing is a 'first order approximation' to a bootstrap, hence why I see this solution as suboptimal.
