%Time Series. Econ 232C. homework 2 problem 1 Anton Cheremukhin, Christoph Bertsch and Martin Szydlowski
close all;
clear all;
clc;
% get data
[SP] = textread('sp500.txt','%f');
T=size(SP,1);
r=log(SP(2:T))-log(SP(1:T-1));
r2=r.*r;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ACF-PACF for returns
%compute ACF
h=r'-mean(r);    T=size(h,2);    L=24;
gamma0=(h*h')/T;
for j=1:L
    gamma(j)=( h(1,(j+1):T) * (h(1,1:(T-j)))')/(T-j);
    ster(j)=2/sqrt(T);
end
ACF=gamma/gamma0;   
figure(1);  subplot(1,2,1); stem_hadles=stem(1:L,ACF,'fill');
hold on;    plot_hadles=plot(1:L,ster,'r',1:L,-ster,'r');   title('ACF r');   hold off;

%compute PACF
for j=1:L
    Y=h(1,(j+1):T)';        X=zeros((T-j),j);
    for k=1:j
        X(:,k)=h(1,(j+1-k):(T-k))';
    end
    beta=inv(X'*X)*(X'*Y);      alpha(j)=beta(j);       ster(j)=2/sqrt(T);
end
PACF=alpha;     subplot(1,2,2);     stem_hadles=stem(1:L,PACF,'fill');
hold on;    plot_hadles=plot(1:L,ster,'r',1:L,-ster,'r');   title('PACF r');  hold off;

%ACF-PACF for hours
%compute ACF
h=r2'-mean(r2);    T=size(h,2);    L=24;
gamma0=(h*h')/T;
for j=1:L
    gamma(j)=( h(1,(j+1):T) * (h(1,1:(T-j)))')/(T-j);
    ster(j)=2/sqrt(T);
end
ACF=gamma/gamma0;   
figure(2);  subplot(1,2,1); stem_hadles=stem(1:L,ACF,'fill');
hold on;    plot_hadles=plot(1:L,ster,'r',1:L,-ster,'r');   title('ACF r^2');   hold off;

%compute PACF
for j=1:L
    Y=h(1,(j+1):T)';        X=zeros((T-j),j);
    for k=1:j
        X(:,k)=h(1,(j+1-k):(T-k))';
    end
    beta=inv(X'*X)*(X'*Y);      alpha(j)=beta(j);       ster(j)=2/sqrt(T);
end
PACF=alpha;     subplot(1,2,2);     stem_hadles=stem(1:L,PACF,'fill');
hold on;    plot_hadles=plot(1:L,ster,'r',1:L,-ster,'r');   title('PACF r^2');  hold off;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ARCH-LM test 
P=4;
h=r(:,1)'-mean(r); T=size(h,2);    
Y=h(1,(P+1):T)';        X=zeros((T-P),P);   ON=ones(size(Y));
for k=1:P
    X(:,k)=(h(1,(P+1-k):(T-k)).^2)';
end
X=[ON X];
beta=inv(X'*X)*(X'*Y);
V=(((Y-X*beta)'*(Y-X*beta))/(size(X,1)-size(X,2)))*inv(X'*X);
GE=eye(size(X,2));  ste=sqrt(spdiags((GE*V*GE').*eye(size(X,2))));
t=GE*beta./ste; pv=2-2*tcdf(abs(t),size(X,1)-size(X,2));
disp('     c        u(-1)       u(-2)       u(-3)       u(-4)   ');    
disp([beta ste t pv]'); disp('R-squared');  R_2=zeros(1,2);
R_2(1,1)=1-((Y-X*beta)'*(Y-X*beta))/((Y-(Y'*ON)/(ON'*ON))'*(Y-(Y'*ON)/(ON'*ON)));
R_2(1,2)=1-(1-R_2(1,1))*(size(X,1)-1)/(size(X,1)-size(X,2));    disp(R_2(1));
disp(T*R_2(1));
disp(['P-value of hypothesis: No serial correlation in squared resduals  ' num2str(1-chi2cdf(T*R_2(1),P))]);
stat(P)=1-chi2cdf(T*R_2(1),P);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%GARCH(1,1)
%use fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options,P1,P2,...) 
e=r(:,1)'-mean(r);      e2=e.*e;        one=1;
options = optimset('Display','off','LargeScale','off','Jacobian','off','MaxFunEvals',100000000,'TolX',0.0000000000001,'TolFun',0.0000000000001,'MaxIter',1000);
%[theta1(i,:),fval1(i,:)] = fmincon(@normlike,[0.25 0.15 0.75],[-1 0 0; 0 -1 0; 0 0 -1; 0 1 1],[-0.00001 -0.0001 -0.0001 one]',[],[],[],[],[],options,e2);
[theta1,fval1] = fmincon(@normlike,[0.25 0.15 0.75],[-1 0 0; 0 -1 0; 0 0 -1; 0 1 1],[-0.00001 -0.0001 -0.0001 one]',[],[],[],[],[],options,e2)
%[theta2(i,:),fval2(i,:)] = fmincon(@normstud,[0.25 0.15 0.75 10],[-1 0 0 0; 0 -1 0 0; 0 0 -1 0; 0 0 0 -1; 0 0 0 1; 0 1 1 0],[-0.00001 -0.0001 -0.0001 -2.0001 100 one]',[],[],[],[],[],options,e2);
[theta2,fval2] = fmincon(@normstud,[0.25 0.15 0.75 10],[-1 0 0 0; 0 -1 0 0; 0 0 -1 0; 0 0 0 -1; 0 0 0 1; 0 1 1 0],[-0.00001 -0.0001 -0.0001 -2.0001 100 one]',[],[],[],[],[],options,e2)
omega=theta1(1);    alpha=theta1(2);    beta=theta1(3);     T=max(size(e2));
s2=zeros(size(e2));     s2(1)=e2(1);
for i=2:T
    s2(i)=omega+alpha*e2(i-1)+beta*s2(i-1);
end
err1=e./sqrt(s2);
omega=theta2(1);    alpha=theta2(2);    beta=theta2(3);     T=max(size(e2));
s2=zeros(size(e2));     s2(1)=e2(1);
for i=2:T
    s2(i)=omega+alpha*e2(i-1)+beta*s2(i-1);
end
err2=e./sqrt(s2);
figure(3)
plot(2:T,err1(2:T),2:T,e(2:T))

%ACF-PACF for error term
%compute ACF
h=err1(2:T)-mean(err1(2:T));    T=size(h,2);    L=24;
gamma0=(h*h')/T;
for j=1:L
    gamma(j)=( h(1,(j+1):T) * (h(1,1:(T-j)))')/(T-j);
    ster(j)=2/sqrt(T);
end
ACF=gamma/gamma0;   
figure(4);  subplot(1,2,1); stem_hadles=stem(1:L,ACF,'fill');
hold on;    plot_hadles=plot(1:L,ster,'r',1:L,-ster,'r');   title('ACF e/s');   hold off;

%compute PACF
for j=1:L
    Y=h(1,(j+1):T)';        X=zeros((T-j),j);
    for k=1:j
        X(:,k)=h(1,(j+1-k):(T-k))';
    end
    beta=inv(X'*X)*(X'*Y);      alpha(j)=beta(j);       ster(j)=2/sqrt(T);
end
PACF=alpha;     subplot(1,2,2);     stem_hadles=stem(1:L,PACF,'fill');
hold on;    plot_hadles=plot(1:L,ster,'r',1:L,-ster,'r');   title('PACF e/s');  hold off;

%ARCH-LM test 
P=4;
%h=err1(2:T)-mean(err1(2:T)); T=size(h,2);    
Y=h(1,(P+1):T)';        X=zeros((T-P),P);   ON=ones(size(Y));
for k=1:P
    X(:,k)=(h(1,(P+1-k):(T-k)).^2)';
end
X=[ON X];
beta=inv(X'*X)*(X'*Y);
V=(((Y-X*beta)'*(Y-X*beta))/(size(X,1)-size(X,2)))*inv(X'*X);
GE=eye(size(X,2));  ste=sqrt(spdiags((GE*V*GE').*eye(size(X,2))));
t=GE*beta./ste; pv=2-2*tcdf(abs(t),size(X,1)-size(X,2));
disp('     c        u(-1)       u(-2)       u(-3)       u(-4)   ');    
disp([beta ste t pv]'); disp('R-squared');  R_2=zeros(1,2);
R_2(1,1)=1-((Y-X*beta)'*(Y-X*beta))/((Y-(Y'*ON)/(ON'*ON))'*(Y-(Y'*ON)/(ON'*ON)));
R_2(1,2)=1-(1-R_2(1,1))*(size(X,1)-1)/(size(X,1)-size(X,2));    disp(R_2(1));
disp(T*R_2(1));
disp(['P-value of hypothesis: No serial correlation in squared resduals  ' num2str(1-chi2cdf(T*R_2(1),P))]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%rolling window
T=max(size(e2));        Tm=1250;
options = optimset('Display','off','LargeScale','off','Jacobian','off','MaxFunEvals',10000,'TolX',0.000000001,'TolFun',0.000000001,'MaxIter',100);
for i=1:T-Tm+1
    e2t=e2(i:Tm+i-1);
    t(i)=i
    if i<2
        theta0=[0.25 0.15 0.75];
    else
        theta0=theta1(i-1,:);
    end
    [theta1(i,:),fval1(i,:)] = fmincon(@normlike,theta0,[-1 0 0; 0 -1 0; 0 0 -1; 0 1 1],[-0.00001 -0.0001 -0.0001 one]',[],[],[],[],[],options,e2t);
end
figure(5)
subplot(2,2,1)
plot(t,theta1(:,1))
subplot(2,2,2)
plot(t,theta1(:,2))
subplot(2,2,3)
plot(t,theta1(:,3))