clear all;
close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = 300;
beta_teo = [3 -1 2.5 0 -1.5]';
p = length(beta_teo)-1;
sigma = 2;
alpha=0.05;
disturbances = 'gaussian';
phi=0.8;
lag=5;
coef_arch=[0.1 0 0 0 .2];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Data Generating Process
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Xobs = rand(N,p);
X = [ones(N,1) Xobs];
switch disturbances
    case 'gaussian'
        % I.i.d. gaussian disturbances
        eps = randn(N,1);
        y = X*beta_teo+sigma*eps;
    case 'autocorrelated'
        % Autocorrelated disturbances
        eps = 0.1*randn(N,1);
        for j=(1+lag):N
            eps(j)=phi*eps(j-lag)+randn(1);
        end
        y = X*beta_teo+eps;
    case 'heteroscedasticity'        
        eps=sqrt(0.8*X(:,1)+X(:,3).^2).*randn(N,1);
        %eps=[1:N]'.*randn(N,1);
        y = X*beta_teo+eps;
    case 'arch'
        order_arch=length(coef_arch)-1;
        a0 = coef_arch(1);
        ai = coef_arch(2:end);
        eps=randn(N+order_arch,1);
        for j=(order_arch+1):(N+order_arch)
            eps(j)= a0;
            for pp=1:order_arch
                eps(j)=eps(j)+ai(pp)*(eps(j-pp))^2;
            end
            eps(j)=sigma*randn(1)*sqrt(eps(j));
        end
        eps=eps(order_arch+1:end);
        y = X*beta_teo+eps;
    case 'unconstancy'
        % Unconstancy of parameters
        eps = randn(N,1);
        n1 = floor(N*0.2);
        n2 = floor(N*0.4);
        y = zeros(N,1);
        y([1:n1 (n2+1):N]) = X([1:n1 (n2+1):N],:)*beta_teo+sigma*eps([1:n1 (n2+1):N]);
        y([(n1+1):n2]) = X([(n1+1):n2],:)*(beta_teo+[0 0 0 4 0]')+sigma*eps([(n1+1):n2]);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Estimation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
reg1 = ols(y,X);
[beta,bint,resi,rint,stats]  = regress(y,X,alpha);
[beta_glm, dev, stats_glm] = glmfit(Xobs,y,'normal','link','identity','constant','on');

%%%
% Regression coefficients beta
%%%
reg1.beta
beta
beta_glm
stats_glm.beta
% Ricorda inv(A) si puo' calcolare
% cone soluzione del sistema A\I
beta_hat = (X'*X\eye(p+1))*(X'*y);
beta_hat1 = (X'*X\eye(p+1))*X'*y;
beta_hat2 = inv(X'*X)*(X'*y);
beta_hat3 = inv(X'*X)*X'*y;
[q r] = qr(X);
beta_hat4 = (r'*r)\eye(p+1)*(X'*y);
beta_hat5 = (X'*X)\(X'*y);

%%%
% Disturbance variance sigma^2
%%%
reg1.sige
stats(4)
(stats_glm.sfit)^2
e = y - X*beta_hat;
k = p+1;
sigma_hat = e'*e/(N-k)

%%%
% Covariance matrix of beta
%%%
stats_glm.covb
covb = sigma_hat*inv(X'*X)

%%%
% Std of regression coefficients std(beta)
%%%
reg1.bstd
stats_glm.se
beta_std = sqrt(diag(covb))

%%%
% t-statistic of regression coefficients beta/std(beta)
%%%
reg1.tstat
stats_glm.t
beta_tstat = beta_hat./beta_std
% p-values of the t-statistic
stats_glm.p
tstat_pvals = 2*(1-tcdf(abs(beta_tstat),N-k))

%%%
% CI of regression coefficients [beta_lower, beta_upper]
%%%
reg1.bint
bint
ttab = tinv(1-alpha/2,N-k);
[beta_hat-ttab*beta_std beta_hat+ttab*beta_std]

%%%
% Residuals
%%%
resid=reg1.resid;
resid_glm = stats_glm.resid;
res = y-X*beta_hat;
[resid resid_glm res]

%%%
% CI of residuals
% E((eps_t)^2) = sigma^2*(1-ht)
% ht = t-th diagonal element of the 
% hat matrix X*inv(X'*X)*X'
% ht = xt'*inv(X'*X)*xt
%%%
tmp = zeros(N,1);
XX = inv(X'*X);
for j=1:N
    c = (X(j,:))';
    tmp(j)=ttab*sqrt(sigma_hat*(1-c'*XX*c));
end
res_lo = res-tmp;
res_hi = res+tmp; 
[rint(:,1) res_lo rint(:,2) res_hi]

%%%
% Coefficient of multiple correlation R
%%
reg1.rsqr
stats(1)
yhat=X*beta_hat;
r2 = (yhat'*yhat-N*(mean(y))^2)/(y'*y-N*(mean(y))^2)

%%%
% Adjusted coefficient of multiple correlation R_bar
%%
reg1.rbar
r2bar = 1 - (res'*res)/(y'*y-N*(mean(y))^2)*(N-1)/(N-k)

%%%
% F-statistic for testing the joint significance
% of the complete set of regressors
%%%
stats(2)
fstat = r2/(1-r2)*(N-k)/(k-1)
stats(3)
fstat_pval = 1-fcdf(fstat,k-1,N-k)

%%%
% Testing linear hypothesis about beta
% R*beta=r
%%%
% 1. Testing the joint significance
%    of the complete set of regressors
R = [zeros(k-1,1) eye(k-1)];
r = zeros(k-1,1);
q = k-1;
fstat = ((R*beta-r)'*inv(R*inv(X'*X)*R')*(R*beta-r)/q)/(res'*res/(N-k))
fstat_pval = 1-fcdf(fstat,q,N-k)
% 2. Testing beta1=0
R = zeros(1,k); R(2)=1;
r = 0;
q = 1;
fstat = ((R*beta-r)'*inv(R*inv(X'*X)*R')*(R*beta-r)/q)/(res'*res/(N-k))
fstat_pval = 1-fcdf(fstat,q,N-k)
% 3. Testing beta4=0
R = zeros(1,k); R(4)=1;
r = 0;
q = 1;
fstat = ((R*beta-r)'*inv(R*inv(X'*X)*R')*(R*beta-r)/q)/(res'*res/(N-k))
fstat_pval = 1-fcdf(fstat,q,N-k)
% 4. Testing beta2=beta5
R = zeros(1,k); R(2)=1; R(5)=-1;
r = 0;
q = 1;
fstat = ((R*beta-r)'*inv(R*inv(X'*X)*R')*(R*beta-r)/q)/(res'*res/(N-k))
fstat_pval = 1-fcdf(fstat,q,N-k)


%%%
% Testing linear hypothesis about beta
% with restricted and unrestricted regressions
%%%
reg_unrest = ols(y,X);
res_unrest=reg_unrest.resid;
su = res_unrest'*res_unrest;
% 1. Testing the joint significance
%    of the complete set of regressors
q=4;
reg_rest = ols(y,X(:,1));
res_rest=reg_rest.resid;
sr = res_rest'*res_rest;
fstat = ((sr-su)/q)/(su/(N-k))
fstat_pval = 1-fcdf(fstat,q,N-k)
% 2. Testing beta1=0
q=1;
reg_rest = ols(y,[X(:,1) X(:,3:end)]);
res_rest=reg_rest.resid;
sr = res_rest'*res_rest;
fstat = ((sr-su)/q)/(su/(N-k))
fstat_pval = 1-fcdf(fstat,q,N-k)
% 3. Testing beta4=0
q=1;
reg_rest = ols(y,[X(:,1:3) X(:,5:end)]);
res_rest=reg_rest.resid;
sr = res_rest'*res_rest;
fstat = ((sr-su)/q)/(su/(N-k))
fstat_pval = 1-fcdf(fstat,q,N-k)
% 4. Testing beta2=beta5
q=1;
reg_rest = ols(y,[X(:,1) X(:,3:4) X(:,2)+X(:,5)]);
res_rest=reg_rest.resid;
sr = res_rest'*res_rest;
fstat = ((sr-su)/q)/(su/(N-k))
fstat_pval = 1-fcdf(fstat,q,N-k)


%%%
% (In-sample) Prediction
%%%
yhat=reg1.yhat;
y_hat=X*beta_hat;
[yhat y_hat]

%%%
% (Out-of-sample) Prediction
%%%
Nout = 20;
Xout = rand(Nout,p);
Xout = [ones(Nout,1) Xout];
yout = Xout*beta_hat;
tmp = zeros(Nout,1);
XX = inv(X'*X);
for j=1:Nout
    c = (Xout(j,:))';
    tmp(j)=ttab*sqrt(sigma_hat*c'*XX*c);
end
yout_lo = yout-tmp;
yout_hi = yout+tmp;
figure(1)
[x idx]=sort(Xout(:,3));
plot(x,yout(idx),'x-b',x,yout_lo(idx),':b',x,yout_hi(idx),':b');

%%%
% Graphics and Tables
%%%
figure(2)
plt(reg1)
prt(reg1)
figure(3)
rdiag(y,X)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Test of parameter constancy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%
% Recursive estimation
%%%
nk=10;
betas=zeros(N-nk*k,k);
beta_lo=betas;
beta_hi=betas;
sigmas=sigma_hat*ones(N-nk*k,1);
for j=nk*k+1:N
    [betas(j-nk*k,:),bint,re,reint,stats]=regress(y(1:j),X(1:j,:));
    beta_lo(j-nk*k,:)=bint(:,1)';
    beta_hi(j-nk*k,:)=bint(:,2)';
    sigmas(j-nk*k)=stats(4);
end
figure(4);
subplot(2,2,1)
plot(nk*k+1:N,betas(:,1),'b',...
    nk*k+1:N,beta_lo(:,1),'r',nk*k+1:N,beta_hi(:,1),'r')
axis tight
set(gca,'xtick',round(linspace(nk*k+1,N,4)))
subplot(2,2,2)
plot(nk*k+1:N,betas(:,2),'b',...
    nk*k+1:N,beta_lo(:,2),'r',nk*k+1:N,beta_hi(:,2),'r')
axis tight
set(gca,'xtick',round(linspace(nk*k+1,N,4)))
subplot(2,3,4)
plot(nk*k+1:N,betas(:,3),'b',...
    nk*k+1:N,beta_lo(:,3),'r',nk*k+1:N,beta_hi(:,3),'r')
axis tight
set(gca,'xtick',round(linspace(nk*k+1,N,4)))
subplot(2,3,5)
plot(nk*k+1:N,betas(:,4),'b',...
    nk*k+1:N,beta_lo(:,4),'r',nk*k+1:N,beta_hi(:,4),'r')
axis tight
set(gca,'xtick',round(linspace(nk*k+1,N,4)))
subplot(2,3,6)
plot(nk*k+1:N,betas(:,5),'b',...
    nk*k+1:N,beta_lo(:,5),'r',nk*k+1:N,beta_hi(:,5),'r')
axis tight
set(gca,'xtick',round(linspace(nk*k+1,N,4)))

%%%
% One-step ahead prediction errors
%%%
nk=10;
rres = zeros(N-nk*k,1);
rresint = rres;
wt=rres;
for j=nk*k+1:N
    [b,bint,re,reint,stats]=regress(y(1:j-1),X(1:j-1,:));
    xj = (X(j,:))';
    rres(j-nk*k)=y(j)-xj'*b;
    XX = inv(X(1:j-1,:)'*X(1:j-1,:));
    rresint(j-nk*k)=ttab*sqrt(stats(4)*(1+xj'*XX*xj));
    wt(j-nk*k)=(y(j)-xj'*b)/sqrt(1+xj'*XX*xj);
end
figure(5)
plot(nk*k+1:N,rres,'b',nk*k+1:N,-rresint,'r',nk*k+1:N,rresint,'r')
pct = sum(rres<-rresint | rres>rresint)/numel(rres)

%%%
% CUSUM and CUSUMSQ tests of Brown,Durbin and Evans
%%%
cs = cusums(y,X);
figure(6)
wt = cs.rres;
wt2 = recresid(y,X);
wt(1:k)=0;
Wt=cumsum(wt);
Wt=Wt/sqrt(sigma_hat);
plot(1:N,Wt,'-r')
line([k N],[ 0.948*sqrt(N-k)  3*0.948*sqrt(N-k)],'Color','b','LineStyle',':')
line([k N],[-0.948*sqrt(N-k) -3*0.948*sqrt(N-k)],'Color','b','LineStyle',':')
line([k N],[0 0],'Color','k','LineStyle',':')
title('CUSUM with 95% intervals');
legend('cusum','upper95','lower95','Location','NorthWest');
xlabel('Observations');

c0tab=[.475 .50855 .46702 .44641 .42174 .40045 .38294 .36697 .35277 .34022 .32894 .31869 ...
    .30935 .30081 .29296 .28570 .27897 .27270 .26685 .26137 .25622 .25136 .24679 .24245...
    .23835 .23445 .23074 .22721 .2238 .2206 .2175 .2146 .2117 .2090 .2064 .2039 .2014...
    .1991 .1968 .1946];
[coef,r,j]=nlinfit(10:40,c0tab(10:40),@c0fun,[2 2 .5]);

figure(7)
plot(1:N,cs.cusums,'-r')
c0=c0fun(coef,N);
line([k N],[c0 1+c0],'Color','b','LineStyle',':')
line([k N],[-c0 1-c0],'Color','b','LineStyle',':')
line([k N],[0 1],'Color','k','LineStyle',':')
title('CUSUMSQ with 95% intervals');
legend('cusumsq','upper95','lower95','mean value','Location','NorthWest');
xlabel('Observations');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Tests for heteroscedasticity
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%
% Breusch-Pagan Test
%%%
bpagan(y,X);
bp=bpagan(y,X);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Tests for autocorrelated disturbances
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%
% Durbin-Watson test
%%%
reg1.dw
dw = sum((diff(res)).^2)/(res'*res)
% If dw < dw_lo, reject the hypothesis of non autocorrelated residuals
% in favor of of the hypothesis of positive first-order autocorrelation
% If dw > dw_up, do not reject the null hypothesis
% If dw_lo < dw < dw_up, the test is inconclusive
% N = 100, k = 5, dw_lo=1.571, dw_up=1.780

%%%
% Box-Pierce-Liung Statistic
%%%
p = [1:8];
[q_stat, q_pvals] = qstat2(reg1.resid,p);
figure(8)
plot(p,q_pvals,'.')
axis([p(1) p(end) 0 1])
set(gca,'Ytick',[0 0.05 .1:.2:1])
line([p(1) p(end)], [0.05 0.05],'LineStyle',':','Color','r')
xlabel('Lag p'): ylabel('p-values');
title('Box-Pierce-Liung Statistic')

%%%
% Testing for ARCH
%%% 
[arch_stat,arch_pvals]=arch(reg1.resid,p);
figure(9)
plot(p,arch_pvals,'.')
axis([p(1) p(end) 0 1])
set(gca,'Ytick',[0 0.05 .1:.2:1])
line([p(1) p(end)], [0.05 0.05],'LineStyle',':','Color','r')
xlabel('ARCH(p)'): ylabel('p-values');
title('Testing for ARCH')