Applying correction to a time series in Matlab
There is really nothing wrong with using a loop in MATLAB. For decades, they have been teaching us to vectorize our code. But in recent years the differences between vectorized code and the equivalent loop code has been shrinking dramatically. And I've run into examples where vectorized code is actually slower!

Nonetheless, there is a simple improvement to be made to your code. Note that what you implemented is an algorithm of \$O(n^2)\$. This simplified bit:


for i = 1:numel(ts)
   ts(i:end) = ts(i:end) + correction(2);
end



does \$n(n+1)/2\$ additions (hence \$O(n^2)\$). Also the indexing is expensive, since you're copying half the array (on average) every loop iteration.

I suggest that, instead of adding correction(2) to all subsequent elements, you store the "current correction", update it every loop iteration, and add it to the current value only.

Below is the test function I wrote. method1 is your code, method2 is my suggestion. At the top is a function that exercises these two methods, compares their output, and times them:

function test_methods
ts = [1 2 3 4 3 2 1 0 -1 -2 -3];
ts1 = method1(ts);
ts2 = method2(ts);
if any(abs(ts1-ts2)&gt;1e-6)
   error('the methods differ');
end
timeit(@()method1(ts))
timeit(@()method2(ts))

function ts = method1(ts)
threshold = [-3.5 3];
correction = [1 -1];
for i = 1:numel(ts)
   if ts(i) &gt; threshold(2)
      ts(i:end) = ts(i:end) + correction(2);
   elseif ts(i) &lt; threshold(1)
      ts(i:end) = ts(i:end) + correction(1);
   end
end

function ts = method2(ts)
threshold = [-3.5 3];
correction = [1 -1];
current = 0;
for i = 1:numel(ts)
   value = ts(i) + current;
   if value &gt; threshold(2)
      current = current + correction(2);
      value = value + correction(2);
   elseif value &lt; threshold(1)
      current = current + correction(1);
      value = value + correction(1);
   end
   ts(i) = value;
end


For the short example input, these functions both run too fast for accurate timing. I see 3.1 &mu;s and 1.4 &mu;s, my version is only twice as fast as yours. But for larger inputs the differences become more important (I figured that the cumulative sum of a random process would imitate appropriately your drifting variable):

ts = cumsum(randn(1,1000));


0.29 ms and 8.33 &mu;s, an order of magnitude difference.

ts = cumsum(randn(1,100000));


Now I see 1.66 s and 0.812 ms, 3 orders of magnitude difference.

Because method1 is quadratic in the input length, and method2 is linear, the time difference grows quadratically.