%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% " The Shale Revolution and Shifting Crude Dymamics"
% Demand shock loading on crude oil price movement
% Author: Liuren Wu, liurenwu@gmail.com
% Final update: 2019/10/06
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all; format compact;

indm=3;k=3; %Define the skew based on floating vol at k=3% and k=-3%
%outputfolder=['../JAER1/'];
yrs=[2004:2:2017]';
dvs=datenum(yrs,1,1);

%load in preprocessed floating IV data for WTI and SPX options
S=load('../spxoptions/spxfloatingskews.mat','ud','xmh','mh','yhh'); %'nkv','matsv','-mat');
sud=S.ud;
inds=find(sud>=datenum(2004,1,1)&sud<=datenum(2016,4,30));
O=load('../futuresoptions/SKews.mat','ud','xmh','mh','yhh'); %,'nkv','matsv','-mat');
oud=O.ud;indo=find(oud>=datenum(2002,1,1)&oud<=datenum(2016,4,30));
[dd,inda,indb]=intersect(sud(inds),oud(indo));
inds=inds(inda); indo=indo(indb);
xmh=S.xmh;

%Define skew based on \mp 3% percentage strike movement
indp=find(round(100*xmh(:))==-k);
inda=find(xmh(:)==0);
indc=find(round(100*xmh(:))==k);
OSKP5=100*(O.yhh(indo,inda)-O.yhh(indo,indp))/k;
OSKC5=100*(O.yhh(indo,indc)-O.yhh(indo,inda))/k;
OAV=sqrt(O.yhh(indo,inda));
SSKP5=100*(S.yhh(inds,inda)-S.yhh(inds,indp))/k;
SSKC5=100*(S.yhh(inds,indc)-S.yhh(inds,inda))/k;
SAV=sqrt(S.yhh(inds,inda));
OA=OAV(:,indm); SA=SAV(:,indm);
OSK=nanmean([OSKP5(:,indm),OSKC5(:,indm)],2);
SSK=nanmean([SSKP5(:,indm),SSKC5(:,indm)],2);


figure(1);clf;
plot(dd,100*OA,'-',dd,100*SA,'--','LineWidth',4);
set(gca, 'XTick',dvs);
set(gca, 'XTickLabel',cellstr(num2str(yrs)))
axis([dd(1)-100,dd(end)+100,0,100])    
legend('WTI','SPX','Location','NorthWest')
ylabel('ATM implied volatility, %','FontSize',16)
set(gca,'Box','on','LineWidth',2,'FontSize', 20)
%print('-depsc','-r70', [outputfolder,'figSPXCLATMV3mts.eps'])

figure(2);clf;
plot(dd,OSK,'-',dd,SSK,'--',dd,dd*0,'-.','LineWidth',4)
set(gca, 'XTick',dvs);
set(gca, 'XTickLabel',cellstr(num2str(yrs)))
axis([dd(1)-100,dd(end)+100,-.8,.4])    
ylabel('Implied variance skew','FontSize',16)
legend('WTI','SPX','Location','NorthWest')
set(gca,'Box','on','LineWidth',2,'FontSize', 20)
%print('-depsc','-r70', [outputfolder,'figSPXCLSK3pt.eps'])

x=[OA,SA];
[mean(x)
min(x)
max(x)]
ind=(dd>datenum(2011,1,0));
[corr(diff(x(~ind,:)))
corr(diff(x(ind,:)))]
[corr(diff(log(x(~ind,:))))
corr(diff(log(x(ind,:))))]

x=[OSK,SSK];
[mean(x)
min(x)
max(x)]
ind=(dd>datenum(2011,1,0));
[corr(diff(x(~ind,:)))
corr(diff(x(ind,:)))]

TOSK=min(0,OSK);
etad3=(OSK.*OA.^2)./(SSK.*SA.^2);
etad=sign(etad3).*(abs(etad3).^(1/3));
etad3=(TOSK.*OA.^2)./(SSK.*SA.^2);
etadt=((etad3).^(1/3));

RC2=(etad.^2).*SA.^2./(OA.^2);
RC=((TOSK.*SA./(SSK.*OA)).^2).^(1/3);
Tu=length(dd);
LW=21;meta=etad*NaN;mrc=meta;mht=meta;metat=meta;
 for t=LW:Tu
     meta(t)=nanmean(etad(t-LW+1:t));metat(t)=nanmean(etadt(t-LW+1:t));
     mrc(t)=nanmean(RC(t-LW+1:t));
 end
 
Ad=SA.^2;
Ao=OA.^2;
Sd=SSK; So=OSK;
T=size(OA,1);
B=[ones(T,1), Ad]\Ao;
C=[ones(T,1), Ad.*Sd./Ao]\So;

ed=nanmean([sqrt(B(2)), C(2).^(1/3)]);
RC3=ed.^2.*Ad./Ao;

figure(3);clf;
plot(dd,meta,'-',dd,dd*0,'-.','LineWidth',4)
set(gca, 'XTick',dvs);
set(gca, 'XTickLabel',cellstr(num2str(yrs)))
axis([dd(1)-100,dd(end)+100,-2,3])    
ylabel('Loading on demand shocks','FontSize',16)
set(gca,'Box','on','LineWidth',2,'FontSize', 20)
%print('-depsc','-r70', [outputfolder,'figSPXCLetad.eps'])

figure(4);clf;
plot(dd,100*mrc,'-','LineWidth',4)
set(gca, 'XTick',dvs);
set(gca, 'XTickLabel',cellstr(num2str(yrs)))
axis([dd(1)-100,dd(end)+100,0,100])    
ylabel('Variance contribution from demand shocks, %','FontSize',16)
set(gca,'Box','on','LineWidth',2,'FontSize', 20)
%print('-depsc','-r70', [outputfolder,'figSPXCLRC.eps'])


