function y = semistat(npoints, alpha, drift)
%y = semistat(npoints, alpha, beta, drift)
% generate a semistationary AR(1) singal.
% 
% npoints - number of sample points
% alpha - noise amplitude in stationary part I(0.95)
% beta - noise amplitude in unstationary part I(1)
% drift - contribution of unstationary part 
%           y=stationary+drift*unstationary

rho_s = 0.95; % stationary
rho_d = 1; %drift
s = zeros(npoints, 1);
d = zeros(npoints, 1);

err_s = randn(npoints,1);
err_d = randn(npoints,1);

for i = 2:npoints
  s(i) =  rho_s*s(i-1)+alpha*err_s(i); 
  d(i) =  rho_d*d(i-1)+alpha*err_d(i);  
end

y = (1-drift)*s+drift*d;