clc 
close all
clear all

load('dataSpeed')

%% ARIMA parameters 
n = 2;          % n initial parameter 
tmp_n = 1;      
e1 = 0;         
e2 = 0;
e3 = 0;
theta1 = 0.30;  % moving average parameters, it can take others values 
theta2 = 0.35;
theta3 = 0.35;
sigma = 0; % YOU HAVE TO TUNE THIS PARAMETER!
uplimit = 200;  
downlimit = 50;


estimated_speed = [];
incidens = [];


%% ARIMA algoritmus

for i=1:length(Avg_speed_per_interval)
   if i < 5
        if (isnan(Avg_speed_per_interval(i))) % if there is no vehicle in the tested interval?
           estimated_speed(i+1) = 0; 
        else
           current_speed = Avg_speed_per_interval(i);
           estimated_speed(i) = current_speed;
           e3=e2;
           e2=e1;
           e1=estimated_speed(i) - current_speed; 
           incidens(i) = 0;
           estimated_speed(i+1) = current_speed; 
        end
           
   else
       if (isnan(Avg_speed_per_interval(i))) % if there was cars in the interval
           incidens(i) = 0; 
           
           estimated_speed(i+1) = estimated_speed(i);
       else
           current_speed = Avg_speed_per_interval(i);                
           if (current_speed > uplimit || current_speed < downlimit) 
                                            
              incidens(i) = 1;
              
           else
               incidens(i) = 0;
               
           end
            
           % következ? mérési pontra becslés
           estimated_speed(i+1) = current_speed - e1*theta1 - e2*theta2 -e3*theta3;
           
           
           tmp_n = Number_of_cars_per_interval(i)/current_speed; 
           n = 1/tmp_n
           
           
           uplimit = estimated_speed(i+1) + n*sigma;
           downlimit = estimated_speed(i+1) - n*sigma;

           e3=e2;
           e2=e1;
           e1=estimated_speed(i)-current_speed ; 
       end
   end
end

%% Number of found incidents
count_ones = sum(incidens(:) == 1);
disp(['Number of incidents: ', num2str(count_ones)]);

%% Incident plo
figure;
subplot(3,1,1);
plot(datetime(time(2:length(time)),'ConvertFrom','datenum'),Avg_speed_per_interval);
xlabel('Time');
ylabel('Speed');
title('Average speed over time intervals');
subplot(3,1,2);
plot(datetime(time(2:length(time)),'ConvertFrom','datenum'),incidens);
xlabel('time');
ylabel('Alarm');
title('Incident occurrence');
subplot(3,1,3);
plot(datetime(time(2:length(time)),'ConvertFrom','datenum'),Number_of_cars_per_interval);
xlabel('time');
ylabel('# of car');
title('Number of vehicles measured in time intervals');