%Version 2, March, 2009.
%These routines can be used by and distributed for non-profit academic purposes
%without any royalty except that the users must cite:
%Perron and Yabu (2009): "Estimating Deterministic Trends with an Integrated or
%Stationary Noise Component,"  Journal of Econometrics  151 (2009), 56-69.
%for any other commercial use and/or comments,  please contact 
%Pierre Perron at perron at bu.edu.
%Even though we tried to make this program error-free we cannot be held
%responsible for any consequences that could result from remaining errors.
%Copyright, Pierre Perron and Tomoyoshi Yabu (2009).
%Traslated from the original version of Gauss to Matlab by Lola Gadea, 2017

load LUS.dat; %read your data (log of US GDP in this case)

dep=LUS(:,1);  %set up the data, dep is the dependent variable

T=length(dep);         %number of sample

kmax=fix(12*(T/100)^(1/4)); %maximum number of lag length

a_type = 2; %type=1; supper-efficient estimator is based on the median unbiased estimator
          %type=2; supper-efficient estimator is based on the upper-biased estimator

%***********************************/
%*Do not change the following code */
%***********************************/
[tstat1,stdb1,b1,ams,tstat2,stdb2,b2,aus,ahat,am,au,p]=tcode(dep,kmax);

fprintf('************************************************\n');
fprintf('The model considered is the following:\n');
fprintf('y[t] = b0 + b1*t + u[t]\n');
fprintf('where u[t] = a*u[t-1] + A(L)(u[t-1]-u[t-2]) + e[t]\n');

if a_type==1;
fprintf('\n');
fprintf(strcat('y[t] =', num2str(b1(1),'%6.4f'),'+',num2str(b1(2),'%6.4f'),'*t + u(t)\n'));
fprintf(strcat(32,32,32,32, '(',num2str(stdb1(1),'%6.4f'),')',32,'(',num2str(stdb1(2),'%6.4f'),')','\n'));
fprintf ('\n');
fprintf('where standard errors are given in parentheses and\n'); 
fprintf (strcat('lag length is chosen as k =',num2str(p,'%6.0f'), 32, 'by Modified AIC\n'));
fprintf ('\n');
fprintf ('weighted symmetric least-squares esimtator: %6.4f\n',ahat)
fprintf ('median-unbiased esimtator(a_MU); %6.4f\n',am) 
fprintf ('supper-efficient estimate(based on a_MU): %6.4f\n',ams)
fprintf ('\n');
fprintf ('t-stat: %6.4f\n',tstat1); 
fprintf ('\n');
display('90% confidence interval for b1:')
fprintf (strcat(num2str(b1(2),'%6.4f'),'-1.65*',num2str(stdb1(2),'%6.4f'), '<',  num2str(b1(2),'%6.4f'),  '<', num2str(b1(2),'%6.4f'),'+1.65*',num2str(stdb1(2),'%6.4f'),'\n'));
display('95% confidence interval for b1:')
fprintf (strcat(num2str(b1(2),'%6.4f'),'-1.96*',num2str(stdb1(2),'%6.4f'), '<',  num2str(b1(2),'%6.4f'),  '<', num2str(b1(2),'%6.4f'),'+1.96*',num2str(stdb1(2),'%6.4f'),'\n'));
display('99% confidence interval for b1:')
fprintf (strcat(num2str(b1(2),'%6.4f'),'-2.57*',num2str(stdb1(2),'%6.4f'), '<',  num2str(b1(2),'%6.4f'),  '<', num2str(b1(2),'%6.4f'),'+2.57*',num2str(stdb1(2),'%6.4f'),'\n'));


elseif a_type==2
fprintf('\n');
fprintf(strcat('y[t] =', num2str(b2(1),'%6.4f'),'+',num2str(b2(2),'%6.4f'),'*t + u(t)\n'));
fprintf(strcat(32,32,32,32, '(',num2str(stdb2(1),'%6.4f'),')',32,'(',num2str(stdb2(2),'%6.4f'),')','\n'));
fprintf ('\n');
fprintf('where standard errors are given in parentheses and\n'); 
fprintf (strcat('lag length is chosen as k =',num2str(p,'%6.0f'), 32, 'by Modified AIC\n'));
fprintf ('\n');
fprintf ('weighted symmetric least-squares estimator: %6.4f\n',ahat)
fprintf ('median-unbiased esimtator(a_UB); %6.4f\n',au) 
fprintf ('supper-efficient estimate(based on a_UB): %6.4f\n',aus)
fprintf ('\n');
fprintf ('t-stat: %6.4f\n',tstat2); 
fprintf ('\n');
display('90% confidence interval for b2:')
fprintf (strcat(num2str(b2(2),'%6.4f'),'-1.65*',num2str(stdb2(2),'%6.4f'), '<',  num2str(b2(2),'%6.4f'),  '<', num2str(b2(2),'%6.4f'),'+1.65*',num2str(stdb2(2),'%6.4f'),'\n'));
display('95% confidence interval for b2:')
fprintf (strcat(num2str(b2(2),'%6.4f'),'-1.96*',num2str(stdb2(2),'%6.4f'), '<',  num2str(b2(2),'%6.4f'),  '<', num2str(b2(2),'%6.4f'),'+1.96*',num2str(stdb2(2),'%6.4f'),'\n'));
display('99% confidence interval for b2:')
fprintf (strcat(num2str(b2(2),'%6.4f'),'-2.57*',num2str(stdb2(2),'%6.4f'), '<',  num2str(b2(2),'%6.4f'),  '<', num2str(b2(2),'%6.4f'),'+2.57*',num2str(stdb2(2),'%6.4f'),'\n'));
end
 


