%	PROGRAMA QUE PERMITE CALCULAR EL TEST DE MARONNA YOHAI 
%	LA SERIE DE REFERENCIA ES LA DE LA VARIABLE X Y LA SERIE 
%       REAL EN LA QUE SE
%	QUIERE IDENTIFICAR LA OCURRENCIA DE UN CAMBIO ES LA DE 
%       LA VARIABLE Y.

clear

Mar=xlsread('M.xls');
[N,L] = size(Mar);

PX=Mar(:,2);
PY=Mar(:,3);
prom1=mean(PX);
prom2=mean(PY);
des1=std(PX);
des2=std(PY);

%   estandarizo las series 

X=(PX-prom1)/des1;
Y=(PY-prom2)/des2;

%	SE GENERA LA SERIE DE VALORES ACUMULADOS
	XX(1)=X(1);
	YY(1)=Y(1);
	XAC(1)=X(1);
	YAC(1)=Y(1);

for J=2:N;
	YAC(J)=YAC(J-1)+Y(J);
	XAC(J)=XAC(J-1)+X(J);
   	XX(J)=XAC(J)/J;
	YY(J)=YAC(J)/J;
end
%	SE CALCULAN LOS VALORES MEDIOS DE AMBAS VARIABLES
	XXM=XX(N);
	YYM=YY(N);
%	SE CALCULAN LAS VARIANZAS Y LA COVARIANZA

SX=var(X)*(N-1);
SY=var(Y)*(N-1);
SXYA=cov(X,Y)*(N-1);
SXY=SXYA(1,2);

%	SE CALCULAN LOS PARÁMETROS QUE PERMITEN IDENTIFICAR EL AÑO
%	POSTERIOR AL CAMBIO EN LA MEDIA (T(L)) Y LA MAGNITUD DEL CAMBIO (D(L))

for L=1:N-1;
	F(L)=SX-(((XX(L)-XXM)^2)*N*L)/(N-L);
	D(L)=(SX*(YYM-YY(L))-SXY*(XXM-XX(L)))*N/((N-L)*F(L));
	T(L)=(L*((N-L)*D(L)^ 2)*F(L))/(SX*SY-SXY^2);
end

plot(Mar(1:N-1,1),T,'r');
xlabel('años');
ylabel('T');
