function [out1,out2,out3,out4]=regress_tvp(Y,X,T0,D0in,V0in,scalex,TT0,reps,burn,varargin)

if isempty(varargin)
    remove_train=1;
else
    remove_train=varargin{1};
end

N=cols(X);


%step 1 set starting values and priors using a pre-sample of 10 quarters
y0=Y(1:T0,:);
x0=X(1:T0,:);
b0=x0\y0;
e0=y0-x0*b0;
sigma0=(e0'*e0)/T0;
V0=sigma0*inv(x0'*x0); 

%prior for the variance of the observation equation p(R)~IG(T0/2, D0/2)
D0=D0in;
v0=V0in;
%priors for the variance of the transition equation
Q0=V0*T0*scalex;  %prior for the variance of the transition equation error
P00=V0;                    % variance of the intial state vector  variance of state variable p[t-1/t-1]
beta0=b0';       % intial state vector   %state variable  b[t-1/t-1]

%initialise 
Q=Q0;
R=sigma0;


%remove intial Sample
if remove_train==1
Y=Y(T0+1:end,:);
X=X(T0+1:end,:);
end
T=rows(X);

%Gibbs sampling algorithm Step 2

mm=1;
for m=1:reps

%%Step 2a Set up matrices for the Kalman Filter

ns=cols(beta0);
F=eye(ns);
mu=0;
beta_tt=[];          %will hold the filtered state variable
ptt=zeros(T,ns,ns);    % will hold its variance



% %%%%%%%%%%%Step 2b run Kalman Filter
i=1;
x=X(i,:);
%Prediction
beta10=mu+beta0*F';
p10=F*P00*F'+Q;
yhat=(x*(beta10)')';                                                
eta=Y(i,:)-yhat;
feta=(x*p10*x')+R;
%updating
K=(p10*x')*inv(feta);
beta11=(beta10'+K*eta')';
p11=p10-K*(x*p10);
beta_tt=[beta_tt;beta11];
ptt(i,:,:)=p11;

for i=2:T
   x=X(i,:);


    %Prediction
beta10=mu+beta11*F';
p10=F*p11*F'+Q;
yhat=(x*(beta10)')';                                               
eta=Y(i,:)-yhat;
feta=(x*p10*x')+R;
%updating
K=(p10*x')*inv(feta);
beta11=(beta10'+K*eta')';
p11=p10-K*(x*p10);

ptt(i,:,:)=p11;
beta_tt=[beta_tt;beta11];

end


%%%%%%%%%%%end of Kalman Filter%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%step 2c Backward recursion to calculate the mean and variance of the distribution of the state
%vector

beta2 = zeros(T,ns);   %this will hold the draw of the state variable
wa=randn(T,ns);
error=zeros(T,1);

i=T;  %period t
p00=squeeze(ptt(i,:,:)); 
beta2(i,:)=beta_tt(i:i,:)+(wa(i:i,:)*chol(p00));   %draw for beta in period t from N(beta_tt,ptt)
error(i,:)=Y(i,:)-X(i,:)*beta2(i:i,:)';  %var residuals


%periods t-1..to .1


for i=T-1:-1:1
pt=squeeze(ptt(i,:,:));

bm=beta_tt(i:i,:)+(pt*F'*inv(F*pt*F'+Q)*(beta2(i+1:i+1,:)-beta_tt(i,:)*F')')';  %update the filtered beta for information contained in beta[t+1] 
                                                                                %i.e. beta2(i+1:i+1,:) eq 8.16 pp193 in Kim Nelson

pm=pt-pt*F'*inv(F*pt*F'+Q)*F*pt;  %update covariance of beta

beta2(i:i,:)=bm+(wa(i:i,:)*chol(pm));  %draw for beta in period t from N(bm,pm)eq 8.17 pp193 in Kim Nelson
error(i,:)=Y(i,:)-X(i,:)*beta2(i:i,:)';  %var residuals

end

% step 3 sample Q from the IW distribution
errorq=diff(beta2);
scaleQ=(errorq'*errorq)+Q0;
Q=iwishrnd(scaleQ,T+TT0);

%step4 sample R from the IG distribution
T1=v0+T;
D1=D0+error'*error;
%draw from IG
z0=randn(T1,1);
z0z0=z0'*z0;
R=D1/z0z0;
disp(strcat('iter',num2str(m)))
if m>burn
    %save output from Gibbs sampler
    out1(mm,1:T,:)=beta2;
    out2(mm,1:1,1:1)=R;
    out3(mm,1:N,1:N)=Q;
    out4(mm,:)=beta2(:,1).^2./(beta2(:,1).^2+R);
    mm=mm+1;
end

end