% PROGRAM USING RECURSIVE RESIDUALS TO FIND STRUCTURAL BREAK 
clc;

% PARAMETERS 
sig = 5;
b1 = [1;0.2;1];
b2 = [1;0.5;1];
T = 100;
T1 = 50;

% VARIABLES 
y = zeros(T,1);
constant = ones(T,1);
trend = (1:1:(1+1*(T-1)))'; 
x = 5*rand(T,1);
xmat = [constant,trend,x];

% SIMULATE DATA BEFORE THE BREAK 
for tt  = (1:1:T1)
	y(tt) = xmat(tt,:)*b1 + sig*randn(1,1); 
end

% SIMULATE DATA AFTER THE BREAK 
for tt  = (T1+1:1:T)
	y(tt) = xmat(tt,:)*b2 + sig*randn(1,1); 
end

% CALCULATE RECURSIVE RESIDUALS 
rr = zeros(T,1);
w = zeros(T,1);
for tt  = (10:1:T)
	b = inv(xmat(1:tt-1,:)'*xmat(1:tt-1,:))*(xmat(1:tt-1,:)'*y(1:tt-1));
	rr(tt) = y(tt) - xmat(tt,:)*b;
	s2 = (y(1:tt-1)-xmat(1:tt-1,:)*b)'*(y(1:tt-1)-xmat(1:tt-1,:)*b)/(tt-1-3);
	w(tt) = rr(tt)/sqrt(s2*(1+xmat(tt,:)*inv(xmat(1:tt-1,:)'*xmat(1:tt-1,:))*xmat(tt,:)'));
end

% THRESHOLDS 
lb = -2.64*ones(T,1);
ub = 2.64*ones(T,1);

% PLOT RECURSIVE RESIDUALS & THRESHOLDS 
hold on
plot(trend,lb);
plot(trend,ub);
plot(trend,w,'ko');