%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SPX option processing
% generate skew at fixed moneyness k and time to maturity tau
% Author: Liuren Wu, liurenwu@gmail.com
% Final update: 2019/10/06
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;
clear all; format short; format compact

%Load in original data, aggregated into mat files, with pre-defined sequence
data='specify data folder'

load([data,'../data/SPXChain19960104-20090323.mat']);
%chain
%Date  expiry  ID     Days    F  S  R  Div K CP IV  Bid Ask Volume Delta gamma Vega Theta .. OptionIntests
% 1      2      3      4      5  6  7  8   9 10 11   12 13    14   15    16    17    18         24
chain=chain(:,[1:13,14,24]);
chain(chain(:,11)<0,11)=NaN;
chain=chain(chain(:,1)<datenum(2009,3,23),:);
ldate=max(chain(:,1));
load([data,'SPXChain20090323-20110504.mat']);
%Date  expiry  ID     Days    F  S  R  Div K CP IV  Bid Ask Volume OptionIntests
% 1      2      3      4      5  6  7  8   9 10 11   12 13    14   15
A=load([data,'SPXChain20110505-20150130.mat'],'nchain','pos','-mat');
B=load('../spxoptions/SPXChain20150202-20160429.mat','nchain','pos','-mat');
IVS=[chain;nchain;A.nchain;B.nchain];


%filering
ind=IVS(:,4)>7&IVS(:,11)>.05&IVS(:,11)<2;
IVS=IVS(ind,:);
dd=IVS(:,1);
ud=unique(dd);
T=length(ud);

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=20;
for t=1:T
    indt=dd==ud(t);
    IVt=IVS(indt,:);
    matd=IVt(:,4);
    umat=unique(matd);nm=length(umat);
    nmc=NaN(nm,1);yhm=NaN(nm,nxh);
    for m=1:nm
        %m
        indm=matd==umat(m);tau=umat(m)/365;
        Nm=sum(indm);nmc(m)=Nm;
        if Nm>npoints
            IVm=IVt(indm,:);
            x=log(IVm(:,9)./IVm(:,5));
            y=IVm(:,11).^2;
            miv=nanmean(sqrt(y));
            z=x./(miv.*sqrt(tau));
            wpc=normcdf(z);
            pc=IVm(:,10);
            wpc(pc==2)=1-wpc(pc==2);
            yy=lqIVkregress(x,y,x,wpc,bds/2);
            e=y-yy; ve=var(e);
            we=wpc.*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/figSPXFoptionNumberMaturities.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/figSPXFoptionNumberStrikes.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/figSPXFoptionAverageMaturity.eps'])

save('../data/spxfloatingskews.mat','ud','xmh','mh','yhh','nkv','matsv','-mat');

