clear all;
close all;
clc;

%% Setup

%Figure grid
global nrows ncols fig
nrows = 3; % Number of figure rows
ncols = 3; % Number of figure columns
fig = 0; % Initialize figure counter
warning('off','MATLAB:print:FigureTooLargeForPage');

% Figure parameters
lw = 1; % line width
ms = 250; % marker size
fs = 11; % font size

%Set up figure
%hold on;
set(gca,'TickLabelInterpreter','latex');
set(gca,'FontSize',fs);
set(gca,'Clipping','off');
set(gca,'Position',[0.10 0.15 0.80 0.70]);
set(gca,'TickDir','out');

%Import data
RGDPall = importdata('GDPalldata.xls','',1);
GDPall = RGDPall.data;
lnRGDP = log(GDPall);
GDP = GDPall(177:292); %Create vector of GDP values that only include data
                        %from 1991-2019
lnGDP = log(GDP);
% Create vector that aligns the dates to the data
time = (1991.0:0.25:2019.75)';
T = size(time,1);
k = 2;

%% Regression
constant = ones(T,1);
%x = sig*randn(T,1);
xmat = [constant,time];
b = inv(xmat'*xmat)*(xmat'*lnGDP);
e = lnGDP - xmat*b;
%sig = sqrt(e'*e)/(T-k);
sig = 0.03;

%% Monte Carlo Experiment of Recursive Residuals

% Parameters
b1 = [1;0.0075];
b2 = [1;0.0074];
r = (1:T)';

% Variables 
y = zeros(T,1);

% Simulate data before the break
for tt  = (1:1:74)
	y(tt) = xmat(tt,:)*b1 + sig*randn(1,1); 
end

% Simulate data after the break
for tt  = (74+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  = (1: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 and thresholds
newfig;
hold on
plot(r,lb,'k');
plot(r,ub,'k');
plot(r,w,'k-');
%Create axes labels and title
xe = 0.01*(2019.75-1991.0);
ye = 0.01*(0.1+0.1);
text(120+xe,-7,'Quarters',...
    'HorizontalAlignment','left','VerticalAlignment','middle','Interpreter','LaTeX');
text(0,3+ye,'Residuals',...
    'HorizontalAlignment','left','VerticalAlignment','bottom','Interpreter','LaTeX');
fs = 13;
text(60,3+10*ye,...
    'Monte Carlo Experiment of Recursive Residuals of Log of Real GDP from $1991Q1 - 2019Q4$',...
    'HorizontalAlignment','center','VerticalAlignment','bottom',...
    'Interpreter','LaTeX','FontSize',fs);
hold off

%Plot recursive residuals and thresholds against time
newfig;
hold on
plot(time,lb,'k');
plot(time,ub,'k');
plot(time,w,'k-');
%Create axes labels and title
xe = 0.01*(2019.75-1991.0);
ye = 0.01*(0.1+0.1);
text(2020.0+xe,-7,'Time',...
    'HorizontalAlignment','left','VerticalAlignment','middle','Interpreter','LaTeX');
text(1990.0,3+ye,'Residuals',...
    'HorizontalAlignment','left','VerticalAlignment','bottom','Interpreter','LaTeX');
fs = 13;
text(2005,3+10*ye,...
    'Monte Carlo Experiment of Recursive Residuals of Log of Real GDP from $1991Q1 - 2019Q4$',...
    'HorizontalAlignment','center','VerticalAlignment','bottom',...
    'Interpreter','LaTeX','FontSize',fs);
hold off

%% Recursive Residuals

%Calculate the standard deviation of the data: sig
%sig = sqrt(var(lnGDP));
sig = sqrt(e'*e)/(T-k);
%Calculate x, the vector of regressors
%x = sig.*rand(T,1);
%Calculate the X matrix
constant = ones(T,1);
xmat = [constant time];

%Calculate recursive residuals

%Identify number of coefficients, k
k = 2;

%Initialize residuals
e = zeros(T,1);
%Initialize scaled residual
w = zeros(T,1);

for tt = (1:1:T)
    %Calculate coefficient estimates, b
    b = inv(xmat(1:tt-1,:)'*xmat(1:tt-1,:))*(xmat(1:tt-1,:)'*lnGDP(1:tt-1));
    %Calculate residuals, e
    e(tt) = lnGDP(tt)-xmat(tt,:)*b;
    %Calculate sum of squared errors, s2
    s2 = (lnGDP(1:tt-1)-xmat(1:tt-1,:)*b)'*(lnGDP(1:tt-1)-xmat(1:tt-1,:)*b)/(tt-1-k);
    %Calculate the scaled residual, w
    w(tt) = e(tt)/sqrt(1+xmat(tt,:)*inv(xmat(1:tt-1,:)'*xmat(1:tt-1,:))*xmat(tt,:)');
end

%Calculate the error bounds
%Calculation of lower bound, lb
lb = (-2.*sqrt(s2).*sqrt(1+xmat(tt,:)*inv(xmat(1:tt-1,:)'*xmat(1:tt-1,:))*xmat(tt,:)'))...
    *ones(T,1);
%Calculation of upper bound, ub
ub = (2.*sqrt(s2).*sqrt(1+xmat(tt,:)*inv(xmat(1:tt-1,:)'*xmat(1:tt-1,:))*xmat(tt,:)'))...
    *ones(T,1);

%Plot recursive residuals and error bounds
newfig;
hold on
%Plot lower bound
plot(time,lb,'k');
%Plot upper bound
plot(time,ub,'k');
%Plot recursive residuals
plot(time,w,'k-');
%Create axes labels and title
xe = 0.01*(2019.75-1991.0);
ye = 0.01*(0.1+0.1);
text(2020.0+xe,-0.1,'Time',...
    'HorizontalAlignment','left','VerticalAlignment','middle','Interpreter','LaTeX');
text(1990.0,0.1+ye,'Residuals',...
    'HorizontalAlignment','left','VerticalAlignment','bottom','Interpreter','LaTeX');
fs = 13;
text(2005,0.1+ye,'Recursive Residuals of Log of Real GDP from $1991Q1 - 2019Q4$',...
    'HorizontalAlignment','center','VerticalAlignment','bottom',...
    'Interpreter','LaTeX','FontSize',fs);
hold off

%% Recusive Residuals of data from 1947-2020

% Create vector that aligns the dates to the data
Rtime = (1947.0:0.25:2020.5)';
RT = size(Rtime,1);

%Calculate the standard deviation of the data: sig
%sig = sqrt(var(lnGDP));
sig = sqrt(e'*e)/(T-k);
%Calculate x, the vector of regressors
%x = sig.*rand(RT,1);
%Calculate the X matrix
constant = ones(RT,1);
xmat = [constant Rtime];

%Calculate recursive residuals

%Identify number of coefficients, k
k = 2;

%Initialize residuals
e = zeros(RT,1);
%Initialize scaled residual
w = zeros(RT,1);

for tt = (1:1:RT)
    %Calculate coefficient estimates, b
    b = inv(xmat(1:tt-1,:)'*xmat(1:tt-1,:))*(xmat(1:tt-1,:)'*lnRGDP(1:tt-1));
    %Calculate residuals, e
    e(tt) = lnRGDP(tt)-xmat(tt,:)*b;
    %Calculate sum of squared errors, s2
    s2 = (lnRGDP(1:tt-1)-xmat(1:tt-1,:)*b)'*(lnRGDP(1:tt-1)-xmat(1:tt-1,:)*b)/(tt-1-k);
    %Calculate the scaled residual, w
    w(tt) = e(tt)/sqrt(1+xmat(tt,:)*inv(xmat(1:tt-1,:)'*xmat(1:tt-1,:))*xmat(tt,:)');
end

%Calculate the error bounds
%Calculation of lower bound, lb
lb = (-2.*sqrt(s2).*sqrt(1+xmat(tt,:)*inv(xmat(1:tt-1,:)'*xmat(1:tt-1,:))*xmat(tt,:)'))...
    *ones(RT,1);
%Calculation of upper bound, ub
ub = (2.*sqrt(s2).*sqrt(1+xmat(tt,:)*inv(xmat(1:tt-1,:)'*xmat(1:tt-1,:))*xmat(tt,:)'))...
    *ones(RT,1);

%Plot recursive residuals and error bounds
newfig;
hold on
%Plot lower bound
plot(Rtime,lb,'k');
%Plot upper bound
plot(Rtime,ub,'k');
%Plot recursive residuals
plot(Rtime,w,'k-');
%Create axes labels and title
xe = 0.01*(2019.75-1991.0);
ye = 0.01*(0.1+0.1);
text(2030.0+xe,-0.3,'Time',...
    'HorizontalAlignment','left','VerticalAlignment','middle','Interpreter','LaTeX');
text(1940.0,0.15+ye,'Residuals',...
    'HorizontalAlignment','left','VerticalAlignment','bottom','Interpreter','LaTeX');
fs = 13;
text(1985,0.15+10*ye,'Recursive Residuals of Log of Real GDP from $1947Q1 - 2020Q3$',...
    'HorizontalAlignment','center','VerticalAlignment','bottom',...
    'Interpreter','LaTeX','FontSize',fs);
hold off

%% Functions 
%Setup of figure grid
function figh = newfig
global nrows ncols fig
%Get the computer's screensize as a vector
screensize = get(groot,'ScreenSize'); 
screenwidth = screensize(3);   
screenheight = 0.975*screensize(4);    
figwidth = screenwidth/ncols;  
figheight = screenheight/nrows;  

%Increment the global figure counter  
fig = fig + 1;  
figh = figure(fig);    
set(figh,'OuterPosition', ...    
    [mod(fig-1,ncols)*figwidth, screenheight - ceil(fig/ncols)*figheight, ...      
    figwidth, figheight]);
end

%savefig.m
%Save a figure as pdf, with whitespace around it minimized
function savefig(figname)  

figh = gcf;  
figh.PaperPositionMode = 'auto';  
fig_pos = figh.PaperPosition;  
figh.PaperSize = [fig_pos(3) fig_pos(4)];  
print(figh, figname,'-dpdf');
end