%  Este programa estima una regresión por umbrales. Aplicación:
%  Estimación de la regresión del tipo contado a 3 meses sobre el tipo forward a 3 meses calculado 
%  3 meses antes.
clear
load tipos.prn -ascii;
r3 = tipos(:,1);
r6 = tipos(:,2);

%Construimos el forward a 3 meses para dentro de 3 meses
forward = 100.*(((1+r6/100).^2)./(1+r3/100)-1.0);
%
T = size(r3,1);
lambda = .95;
  [bhat,residuos,sigma2,vcov,estadt,pvalort,R2]=mco_breve( r3(2:T,1) , [ones(size(r3,1)-1,1) r3(1:T-1,1)]);
%   Al calcular los residuos perdemos una observacion. Definimos T como el numero de residuos.
%   Es una unidad menor que el numero de datos de tipos contado o forward.
T = size(residuos,1); 

%  Calculamos la serie de volatilidades del tipo a 3 meses utilizando el modelo EWMA
sigma(1,1)= sigma2;
for i=2:T
    sigma(i,1) = lambda*sigma(i-1,1)+(1-lambda)*residuos(i-1,1)^2;
end
plot(sigma);

%  Construcción de la amtriz de dats a utilizar en las regresiones
%  Tomamos 63 como el numero medio de datos en un periodo de 3 meses
%  Nótese que los datos desde la observacion 64 a la T+1 en r3 se
%  corresponden con los datos en la serie de volatilidad en las posiciones
%  63 hasta la T
%  Hemos de tomar: el tipo de interés a 3 meses en t; el forward a 3 meses
%  calculado hace 3 meses (t-63), y la variable que vayamos a utilizar como
%  determinante del cambio de régimen.
%  Dicha variable podria ser la volatilidad estimada en t. Alternativamente,
%  podria ser el spread entre los tipos a 3 y 6 meses, observado en t.

 z = sigma;
% z = r6-r3;
datos = [r3(64:T,1) forward(1:T-63,1) z(64:T,1)];

%  Descomponemos la muestra en submuestras de alta y baja volatilidad
T =size(datos,1);

SumaRes = [];
Coeficientes = [];
umbral = .0;
maxi = max(datos(:,3));
mini = min(datos(:,3));
num = 20;  %  Número de puntos en la red de volatilidades en la que vamos a buscar
salto = (maxi-mini)/num;
%   Al especificar el indice a continuacion, terminamos antes de llegar al maximo de volatilidad 
%   en la muestra global para asegurarnos de que ambas submuestras tienen suficientes observaciones 
%   para estimar las regresiones. Por la misma razon, puede ser preciso comenzar en j=2.

%  Estimación estándar, con un sólo régimen
[bhat,residuos,sigma2,vcov,estadt,pvalort,R2] = mco_breve(datos(:,1),[ones(T,1) datos(:,2)]);
SCRR = residuos'*residuos;

for j = 1:num-3   
    umbral = mini+j*salto;
    low = [];
    high = [];
    for i=1:T
        if datos(i,3)<umbral
            low = [low; datos(i,:)];
        else
            high = [high; datos(i,:)];
        end
    end
%
Tlow = size(low,1);
Thigh = size(high,1);
[bhatl,residuosl,sigmal2,vcovl,estadtl,pvalortl,R2l] = mco_breve(low(:,1),[ones(Tlow,1) low(:,2)]);
[bhath,residuosh,sigmah2,vcovh,estadth,pvalorth,R2h] = mco_breve(high(:,1),[ones(Thigh,1) high(:,2)]);
estimaciones = [bhatl estadtl bhath estadth];
Coeficientes = [Coeficientes; estimaciones];
SCR = residuosl'*residuosl+residuosh'*residuosh;
SumaRes = [SumaRes; SCR umbral Tlow Thigh];
end
%
plot(SumaRes(:,2),SumaRes(:,1),'-*')
xlabel('Umbrales posibles de cambio de régimen. Son valores de Z')
ylabel('Suma de Cuadrados de Residuos') 
%
[a,b]=find(SumaRes(:,1)==min(SumaRes(:,1)));
display('SCR-mínimo Umbral óptimo  Tamaño muestra "Baja"  Tamaño muestra "Alta"')
SumaRes(a,:)
display('Coeficientes optimos (primera columna) y sus estadísticos t-Student(segunda columna)')
display('Columnas 1 y 2: Regimen BAJO. Columnas 3 y 4: Regimen ALTO')
Coeficientes(2*a-1:2*a,:)
display('Estadistico F')
((SCRR-SumaRes(a,1))/SCRR)*(T-4)/2

