%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% WTI futrues option processing
% generate skew at fixed moneyness k and time to maturity tau
% Author: Liuren Wu, liurenwu@gmail.com
% Final update: 2019/10/06
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all; format compact

%load in original data, compiled into mat format
load('../futuresoptions/IVSC.mat','IVS','udc','flds','-mat'); %2003/7/18-2017/2/8
%flds={'date','matd','maty','F','r','K','P','pc','IV'}; %data sequence in IVS
T=length(ud)
datevec(ud([1,end]))

%filter
indp=IVS(:,9)>.1&IVS(:,9)<1.5&IVS(:,1)>=datenum(2004,1,1);
IVS=IVS(indp,:);

dd=IVS(:,1);
bds=.5;w0=1;xh=[-.2:.01:.2]';nxh=length(xh);graphs=0;
mh=[30,60,91,182,365]'; nmh=length(mh);
yhh=NaN(T,nmh*nxh);
nkv=NaN(T,100);matsv=nkv;
npoints=10;
for t=1:T
    indt=dd==ud(t);
    IVt=IVS(indt,:);
    matd=IVt(:,2);
    umat=unique(matd);nm=length(umat);
    nmc=NaN(nm,1);yhm=NaN(nm,nxh);
    for m=1:nm
        %m
        indm=matd==umat(m);
        Nm=sum(indm);nmc(m)=Nm;
        if Nm>npoints
            IVm=IVt(indm,:);
            x=log(IVm(:,6)./IVm(:,4));
            y=IVm(:,9).^2;
            indc=find(x>0);
            indp=x<=0;
            [xc0,indc0]=min(x(indc));
            yc0=y(indc(indc0));
            yc0h=lqIVkregress(x(indp),y(indp),xc0,w0,bds); %extrapolate put to the first call
            y(indc)=max(.09.^2,y(indc)-yc0+yc0h); %put-call bias correction, use put as benchmark
            yy=lqIVkregress(x,y,x,w0,bds/2);
            e=y-yy; ve=var(e);
            we=exp(-e.^2/(4*ve));
            yh=lqIVkregress(x,y,xh,we,bds);
            if graphs
                figure(1);plot(x,sqrt(y),'o',xh,sqrt(yh),'r','LineWidth',2);grid;xlabel('k','FontSize',16);ylabel('Implied variance','FontSize',16);title(num2str(umat(m)));pause
            end
            yhm(m,:)=yh;
        end
    end
    nkv(t,1:nm)=[nmc];
    matsv(t,1:nm)=[umat];
    inds=nmc>npoints;
    if sum(inds)>=3
        nmc'
        t
        umat=umat(inds);
        yhm=yhm(inds,:);
        xx=interp1([0;umat],[zeros(1,nxh);yhm.*repmat(umat,1,nxh)],mh,'linear','extrap')./repmat(mh,1,nxh);
        yhh(t,:)=xx(:)';
    end
end
xhm=repmat(xh,1,nmh);
xmh=xhm';

figure(1);clf;
plot(ud,nansum(isfinite(nkv),2),'LineWidth',2);
datetick('x','yy');grid
ylabel('Number of Maturities','FontSize',15)
set(gca,'Box','on','LineWidth',2,'FontSize', 16)
%print('-depsc','-r70', ['../paper/figFoptionNumberMaturities.eps'])

figure(2);clf;
plot(ud,nanmean(nkv,2),'LineWidth',2);
datetick('x','yy');grid
ylabel('Number of Strikes Per Maturity','FontSize',15)
set(gca,'Box','on','LineWidth',2,'FontSize', 16)
%print('-depsc','-r70', ['../paper/figFoptionNumberStrikes.eps'])

figure(3);clf;
plot(ud,nanmean(matsv,2)/365,'LineWidth',2);
datetick('x','yy');grid
ylabel('Avrage Option Maturity, Years','FontSize',15)
set(gca,'Box','on','LineWidth',2,'FontSize', 16)
%print('-depsc','-r70', ['../paper/figFoptionAverageMaturity.eps'])

save('../futuresoptions/SKews.mat','ud','xmh','mh','yhh','nkv','matsv','-mat');
return
