Change Detection Using Volatility Filters

gistfile1.txt

#LMS change Detector 
function [state,h]=ChangeDetectionLMSShiftMultipleSensor1(x,WinSlow,WinFast,WinDesired,muBase)
%Slow Response to change, accurate in steady state

[m,n]=size(x);
if m>n
   x=x';
end
[NumSensors,n]=size(x);
shift=300;
xSlow=[zeros(NumSensors,WinSlow),x];
for i=1:length(xSlow)-WinSlow
OutSlow(:,i)=std(xSlow(:,i:i+(WinSlow-1))');
end
OutSlow=[zeros(NumSensors,shift),OutSlow(:,1:end-shift)];
%Fast Response to change, poor accuracy in steady state
xFast=[zeros(NumSensors,WinFast),x];
for i=1:length(xFast)-WinFast
OutFast(:,i)=std(xFast(:,i:i+(WinFast-1))');
end
%Fast Response to change, poor accuracy in steady state
xDesired=[zeros(NumSensors,WinDesired),x];
for i=1:length(xDesired)-WinDesired
OutDesired(:,i)=std(xDesired(:,i:i+(WinDesired-1))');
end
%Adaptive Filter sparse LMS convex combination of two variance estimation
%algorithms
epsilon=0.001;
omega=0;
state=zeros(1,length(x));
i=1;
future=WinSlow+shift+50;
for i=1:NumSensors
    mu(i)=muBase/var(x(1:future)); 
end
count=1;
flag=0;
i=1;

W(1,1:NumSensors)=0;

 while i<length(OutSlow)-WinDesired
        WTemp=mean(W(i,:));
        h(i)=WTemp;

     for ns=1:NumSensors   
        yEstimate=WTemp*OutFast(ns,i) + (1-WTemp)*OutSlow(ns,i);
        %yOut(i)=yEstimate;
                variation=(randn(1)*epsilon);
        
        e(i)=OutDesired(ns,i+WinDesired)-yEstimate;
        %variation=epsilon;
        W(i+1,ns)=WTemp + mu(ns)*e(i)*(abs(WTemp)-variation)*(OutFast(ns,i)-OutSlow(ns,i));
        if W(i+1,ns)>1
            W(i+1,ns)=1;
        end
        if W(i+1,ns)<0
            W(i+1,ns)=0;
        end
           
        
     end
     
     if h(i)>0.8  && flag==0  && i<(length(x)-future)
            if state(i)==1
                state(i+2:i+future)=1;
            else
                state(i:i+future)=1;     
            end
            for f=1:NumSensors
                mu(f)=muBase/var(x(f,i:i+future));
            end
            %flag=1;
            flag=1;
            count=1;
        end
        if  flag==1 
            count=count+1;
           % flag=0;
           if count>future
               flag=0;
           end
        end
        i=i+1;
     
    end
    
    
    
##Statistical Change Detector

    function  [ChangePoint,ActChange]=PeakChangeDetector(x,WinSlow,sigLevel)
    
    xSlow=[zeros(1,WinSlow),x];
    for i=1:length(xSlow)-WinSlow
        OutSlow(i)=std(xSlow(i:i+(WinSlow-1)));
    end
    
    %Differencing of the volatility filter output
    DiffVolTimeSeries=tsmom(OutSlow',WinSlow);
    DiffVolTimeSeries(1:WinSlow)=0;
    
    N=length(DiffVolTimeSeries);
    count=1;
    mu=sqrt(2)*gamma((WinSlow+1)/2)/(gamma((WinSlow)/2)); %Estimate of the first moment of the chi distribution 
    ChangePoint=zeros(1,N);
    h=1;
    ActChange=zeros(1,N);
    %Loop to calculate the changePoints
    while  count<N+1   
        if count == 1
            count=WinSlow*2;
            VarEstimate=var(x(1:count));
            StdEstFilter=2*sqrt(VarEstimate*(WinSlow-(mu^2))/WinSlow);
            count=count+1;
        end
        
        if abs(DiffVolTimeSeries(count))>StdEstFilter*sigLevel  && count+WinSlow<N
            VarEstimate=var(x(count:count+WinSlow));
            StdEstFilter=2*sqrt(VarEstimate*(WinSlow-(mu^2))/WinSlow);
            %store(h)=count;
            ActChange(count)=1;
            h=h+1;
            [~,I] = max(abs(DiffVolTimeSeries(count:count+WinSlow)));
            ChangePoint(count+I-1-WinSlow)=1;
            count=count+WinSlow*2;
        else
            count=count+1;
        end
    end
    
##Generating Random Change Point locations
function [x,timeinstants,sigma,time]=randomPieceWiseGeneratorMC(Num,c)

N=randi([10000,40000]);
x=zeros(Num,N);
xrand=x;
init=1;
flag=0;
count=1;
timeinstants=zeros(1,N);

n=Num;
C=(ones(n,n)*c)-(eye(n)*c)+eye(n);
[E,L] = eig(C);


while flag<1
    initnew=randi([init+1000,init+4000]);
    time(count)=initnew;
    timeinstants(initnew)=1;
    if count==1
        sigma(count)=randi([1,10]);
    else
        r=round(rand(1));
        rLow=(rand(1)*0.6)+0.1;
        rHigh=(rand(1)*3)+1.5;
        sigma(count)=sigma(count-1)*(r*rLow + (1-r)*rHigh);
    end
    x(:,init:initnew)=sigma(count)*E*sqrt(L)*randn(Num,initnew-init+1);
    count=count+1;
    if (initnew+4000)>N
        flag=1;
        timeinstants(time(count-1))=0;
    end
    init=initnew;
end
x(:,initnew:N)=E*sqrt(L)*randn(Num,N-initnew+1)*sigma(count-1);
time=time(1:end-1);


##Script 
%% LMS test script 

    [x,timeinstants,sigma,time]=randomPieceWiseGeneratorMC(4,0.1);
    WinDesired=10;     
     [state,h]=ChangeDetectionLMSShiftMultipleSensor1(x,250,20,10,0.3);    
     plot(x')
     hold
     plot(state*std(x(1,:)))
     
     
     %% Statistical Change Detector 
      
     [x,timeinstants,sigma,time]=randomPieceWiseGeneratorMC(1,0);

      [ChangePoint,ActChange]=PeakChangeDetector(x,200,3);
      plot(x)
      hold
      plot(ActChange*std(x),'r') %Detection  of change point 
      plot(ChangePoint*std(x),'g')%Estimate  of change point location