%This script is to generate Fig1 and Fig2 from Section 2.  
%It iustrates various qpproches to multivariate quantiles i.e. no regression, location case is considered . 
%Terminology follows Halinn and Šiman (2017)

clear;clc;close all
n = 10; %number of simulated data points
Y = [normrnd(0,1,n,1),normrnd(0,1,n,1)]';


%% Figure 1: Direct approach - spatial quantiles 

figure
set(gcf, 'Color', 'w');

v = [-0.5;-0.5];
objfun = @(q) spatq(q,n,Y,v);
[q_sp,fval,exitflag] = fminunc(objfun,[0;0]);
subplot(1,3,1)
hold on
a1 = plot(Y(1,:),Y(2,:),'ok');
set(gca, 'XAxisLocation', 'origin', 'YAxisLocation', 'origin')
xlim([-max(max(abs(Y))),max(max(abs(Y)))]);
ylim([-max(max(abs(Y))),max(max(abs(Y)))]);
a4 = plot(q_sp(1),q_sp(2),'*r','LineWidth',1.5);
for i=1:n
   plot([q_sp(1),Y(1,i)],[q_sp(2),Y(2,i)],'k-','LineWidth',0.5); 
end
arrow([0,0],v);
legend([a1,a4],'simulated data','geometric quantile \bf{Q_{u}}^{(10)}');
title('u=[-0.5,-0.5]')


v = [-0.1;-0.1];
objfun = @(q) spatq(q,n,Y,v);
[q_sp,fval,exitflag] = fminunc(objfun,[0;0]);
subplot(1,3,2)
hold on
a1 = plot(Y(1,:),Y(2,:),'ok');
set(gca, 'XAxisLocation', 'origin', 'YAxisLocation', 'origin')
xlim([-max(max(abs(Y))),max(max(abs(Y)))]);
ylim([-max(max(abs(Y))),max(max(abs(Y)))]);
a4 = plot(q_sp(1),q_sp(2),'*r','LineWidth',1.5);
arrow([0,0],v);
title('u=[-0.1,-0.1]')

v = [0.5;-0.5];
objfun = @(q) spatq(q,n,Y,v);
[q_sp,fval,exitflag] = fminunc(objfun,[0;0]);
subplot(1,3,3)
hold on
a1 = plot(Y(1,:),Y(2,:),'ok');
set(gca, 'XAxisLocation', 'origin', 'YAxisLocation', 'origin')
xlim([-max(max(abs(Y))),max(max(abs(Y)))]);
ylim([-max(max(abs(Y))),max(max(abs(Y)))]);
a4 = plot(q_sp(1),q_sp(2),'*r','LineWidth',1.5);
arrow([0,0],v);
title('u=[0.5,-0.5]')




%% Figure 2: Directional quantile 
tau = 0.20;
u = [1/sqrt(2);1/sqrt(2)];      %direction vector, u1*x2 = u2*x1 + k;
Gu = [-1/sqrt(2);1/sqrt(2)];

%projections of data into new coordinate system defined by derections u and Gu
Yu = u'*Y;
Yu_= Gu'*Y;

q_d = quantile(Yu,tau,3); 

%first get coordinates of data projected on u in Y coordinates
Y2 = (-Gu(1,1)*Yu) / (Gu(2,1)*u(1,1) - Gu(1,1)*u(2,1));
Y1 = (Yu - u(2,1)*Y2) / u(1,1);
%and quantile
Y_q = (-Gu(1,1)*q_d) / (Gu(2,1)*u(1,1) - Gu(1,1)*u(2,1));

figure
set(gcf, 'Color', 'w');
hold on
a1 = plot(Y(1,:),Y(2,:),'ok');
set(gca, 'XAxisLocation', 'origin', 'YAxisLocation', 'origin')
xlim([-max(max(abs(Y))),max(max(abs(Y)))]);
ylim([-max(max(abs(Y))),max(max(abs(Y)))]);
ylimits = ylim;
xlimits = xlim;
a4 = plot([xlimits(1,1),xlimits(1,2)],...
    [(u(1,1)/u(2,1))*xlimits(1,1),(u(1,1)/u(2,1))*xlimits(1,2)]);
a3 = plot(Y1,Y2,'xr','LineWidth',1.5);
a5 = plot(Y_q,Y_q,'*g');
legend([a1,a4,a3,a5],'simulated data','direction','projection of data','directional quantile');
title('Directional quantile');




