function solow_detrending(model,savedata,savepicts);
% function solow_detrending(model,savedata,savepicts);
%
% savedata = 1 : saves the data in a MATLAB file
% savepicts = 1 : saves the figures in color EPS files
%
% Last update: April 27, 2001
%
% Author: Klaus Reiner Schenk-Hoppé 
% Email: klaus@iew.unizh.ch
% Web: www.iew.unizh.ch/home/klaus


N = 4.0*150;  % Periods for convergence to random fixed point 
M = 4.0*50;   % Periods for detrending

timeaxis = (0:0.25:M/4.0);
L(1) = 1.0;
for t = 1:M+N,
  L(t+1) = 1.00 * L(t);   % L(t) = 1 for all times
end;


  A = 0.95;
  B = 0.95;
  phi2 = 0.005;
  s = 0.25;

if model==1,
  alpha = 0.75;
  delta = 0.9;
  rho2 = 0.0005;
end
if model==2,
  alpha = 0.25;
  delta = 0.9;
  rho2 = 0.0005;
end
if model==3,
  alpha = 0.25;
  delta = 0.1;
  rho2 = 0.0005;
end

rho1 = - rho2;
phi1 = - phi2;
xi(1) = 0.0;

Erate = 1.0075;
eta(1) = 0.0;

a(1) = 1.0;
z(1) = 1.0;
k(1) = a(1)*(s*z(1)/(eta(1)+delta))^(1.0/(1.0-alpha));

%% generate the sample paths of the stochastic variables first

rst = rand('state');


for t = 1:N+M,
  eta(t+1) = A * eta(t) + phi1 + (phi2-phi1)*rand;
  Irate(t) = Erate + eta(t+1);
  a(t+1) = Irate(t) * a(t);
end

for t = 1:N+M,
  xi(t+1) = B * xi(t) + rho1 + (rho2-rho1)*rand;
  z(t+1) = 1.0 + xi(t+1);
end


%% simulate the stochastic economy

for t = 1:N+M,
  y(t) = z(t) * k(t)^alpha;
  k(t+1) = ((1-delta) * k(t) + s * y(t))/Irate(t);
  interest(t) = alpha * z(t) * k(t)^(alpha-1);
  wage(t) = (1-alpha) * z(t) * k(t)^alpha;
end;


%% calculate quarterly data from periods N to M+N

Labor =  L(N:N+M);
zAggr =  z(N:N+M);
GDP = a(N:N+M) .* Labor .* y(N:N+M);
OutputPerWorker = GDP./Labor;

%NetInterestRate=interest(N+1:N+M).*(N+1:3:N+M-3)+z(N+2:3:N+M-2)+z(N+3:3:N+M);;((1+i1)*(1+i2)*(1+i3)-1
%GrossInterestRate = (1-delta+interest(N+1:N+M)).^4;

CapitalStock = a(N:N+M) .* Labor .* k(N:N+M);
CapitalPerWorker = CapitalStock./Labor;
%WageRate =  a(N:N+M) .* wage(N:N+M);
%OutputPerEffUnit = GDP./(zAggr.*Labor); % per unit of efficient labor (i.e. adjusted for total hours worked)
TechnState = a(N:N+M);
LogTechnState = log(a(N:N+M));
TechnProg = Irate(N:N+M) - 1.0;
OutputIntensity = y(N:N+M);
CapitalIntensity = k(N:N+M);
LogCapitalIntensity = log(k(N:N+M));
LogOutputIntensity = log(y(N:N+M));
ExpLogCapitalIntensity = sum(LogCapitalIntensity)./(M+1);
ExpCapitalIntensity = sum(CapitalIntensity)./(M+1);
ExpOutputIntensity = sum(OutputIntensity)./(M+1);
ExpLogOutputIntensity = sum(LogOutputIntensity)./(M+1);

LogGDP = log(GDP);
LogCapitalStock = log(CapitalStock);
%LogWageRate = log(WageRate);
%LogOutputPerWorker = log(OutputPerWorker);


fig1 = figure;
hold
plot(timeaxis,zAggr,'-b');
plot([timeaxis(1),timeaxis(length(timeaxis))],[1,1],'-k')
grid on
axis('tight')
title('z(t)');

fig2 = figure;
hold
plot(timeaxis,CapitalIntensity,'-b');
plot([timeaxis(1),timeaxis(length(timeaxis))],[ExpCapitalIntensity,ExpCapitalIntensity],'-k')
grid on
axis('tight')
title('Capital intensity k(t)');

fig3 = figure;
hold
plot(timeaxis,OutputIntensity,'-b');
plot([timeaxis(1),timeaxis(length(timeaxis))],[ExpOutputIntensity,ExpOutputIntensity],'-k')
grid on
axis('tight')
title('Output intensity y(t)');

fig4 = figure;
hold
plot(timeaxis,LogCapitalStock,'-b');
plot(timeaxis,ExpLogCapitalIntensity+LogTechnState,'-k');
grid on
axis('tight')
title('Log(Capital stock K(t)) (blue), actual trend (black)');

fig5 = figure;
hold
plot(timeaxis,1+TechnProg,'-b');
plot([timeaxis(1),timeaxis(length(timeaxis))],[Erate,Erate],'-k')
grid on
axis('tight')
title('Rate of technical progress 1+eta_{t+1}');


fig6 = figure;
hold
plot(timeaxis,LogTechnState,'-b');
%plot([timeaxis(1),timeaxis(length(timeaxis))],[1+C,1+C],'-k')
grid on
axis('tight')
title('Log(a(t))');


%%%%%%%%%%%%%% linear trend
%[lr_LogGDP,lr_dataLogGDP] = lin_regress(timeaxis,LogGDP);
%[lr_LogCapitalStock,lr_dataLogCapitalStock] = lin_regress(timeaxis,LogCapitalStock);
%[lr_LogWageRate,lr_dataLogWageRate] = lin_regress(timeaxis,LogWageRate);
%[lr_LogOutputPerWorker,lr_dataLogOutputPerWorker] = lin_regress(timeaxis,LogOutputPerWorker);
%[lr_LogTechnState,lr_dataLogTechnState] = lin_regress(timeaxis,LogTechnState);
%lr_LogGDP
%lr_LogTechnState
%l_detrended_LogGDP = LogGDP - lr_dataLogGDP;
%l_detrended_LogCapitalStock = LogCapitalStock - lr_dataLogCapitalStock;
%l_detrended_LogWageRate = LogWageRate - lr_dataLogWageRate;
%l_detrended_LogOutputPerWorker = LogOutputPerWorker - lr_dataLogOutputPerWorker;
%l_detrended_LogTechnState = LogTechnState - lr_dataLogTechnState;


%%%%%%%%%%%%%% HP filter

HPLogGDP = hpfilter(LogGDP,1600)*transpose(LogGDP);
HPLogCapitalStock = hpfilter(LogCapitalStock,1600)*transpose(LogCapitalStock);
%HPLogWageRate = hpfilter(LogWageRate,1600)*transpose(LogWageRate);
%HPLogOutputPerWorker = hpfilter(LogOutputPerWorker,1600)*transpose(LogOutputPerWorker);
%HPLogTechnState = hpfilter(LogTechnState,1600)*transpose(LogTechnState);

HPtrendLogGDP = LogGDP-transpose(HPLogGDP);
HPtrendLogCapitalStock = LogCapitalStock-transpose(HPLogCapitalStock);
%HPtrendLogWageRate = LogWageRate-transpose(HPLogWageRate);
%HPtrendLogOutputPerWorker = LogOutputPerWorker-transpose(HPLogOutputPerWorker);
%HPtrendLogTechnState = LogTechnState-transpose(HPLogTechnState);

%HPLogWageRate4 = hpfilter(LogWageRate,4)*transpose(LogWageRate);
%HPLogWageRate40 = hpfilter(LogWageRate,40)*transpose(LogWageRate);
%HPLogWageRate400 = hpfilter(LogWageRate,400)*transpose(LogWageRate);


%%%%%%%%%%%%%% BP filter

BPLogGDP = bpf(transpose(LogGDP),6,32,12);
BPLogCapitalStock = bpf(transpose(LogCapitalStock),6,32,12);
%BPLogWageRate = bpf(transpose(LogWageRate),6,32,12);
%BPLogOutputPerWorker = bpf(transpose(LogOutputPerWorker),6,32,12);
%BPLogTechnState = bpf(transpose(LogTechnState),6,32,12);

BPtrendLogGDP = LogGDP-transpose(BPLogGDP);
BPtrendLogCapitalStock = LogCapitalStock-transpose(BPLogCapitalStock);
%BPtrendLogWageRate = LogWageRate-transpose(BPLogWageRate);
%BPtrendLogOutputPerWorker = LogOutputPerWorker-transpose(BPLogOutputPerWorker);
%BPtrendLogTechnState = LogTechnState-transpose(BPLogTechnState);

if savedata == 1,
  if model==1,
     save solow1data;
  end
    if model==2,
     save solow2data;
  end
  if model==3,
     save solow3data;
  end
end;

%HPtmpLogGDP = hpfilter(LogGDP,200000)*transpose(LogGDP);
%HPtmptrendLogGDP = LogGDP-transpose(HPtmpLogGDP);

fig7 = figure;
hold
%plot(timeaxis,HPLogGDP+transpose(HPtrendLogGDP),'-g')
plot(timeaxis,BPtrendLogGDP,'-g')
plot(timeaxis,HPtrendLogGDP,'-m')
%plot(timeaxis,HPtmptrendLogGDP,'-k')
plot(timeaxis,LogGDP,'-b')
plot(timeaxis,ExpLogOutputIntensity+LogTechnState,'-k');
%plot(timeaxis,LogTechnState-lr_LogTechnState(1)+lr_LogGDP(1),'-k')
axis('tight')
title('Trend of log(Y_t): HP1600 (magenta), BK(6,32,12) (green), actual trend (black). Data (blue) for reference')

fig8 = figure;
hold
%plot(timeaxis,HPLogCapitalStock+transpose(HPtrendLogCapitalStock),'-g')
plot(timeaxis,BPtrendLogCapitalStock,'-g')
plot(timeaxis,HPtrendLogCapitalStock,'-m')
plot(timeaxis,LogCapitalStock,'-b')
plot(timeaxis,ExpLogCapitalIntensity+LogTechnState,'-k')
%plot(timeaxis,LogTechnState-lr_LogTechnState(1)+lr_LogCapitalStock(1),'-k')
axis('tight')
title('Trend of log(K_t): HP1600 (magenta), BK(6,32,12) (green), actual trend (black). Data (blue) for reference')

fig9 = figure;
hold
plot(timeaxis,BPLogGDP,'-g')
plot(timeaxis,HPLogGDP,'-m')
plot(timeaxis,LogGDP-(ExpLogOutputIntensity+LogTechnState),'-b')
plot([timeaxis(1),timeaxis(length(timeaxis))],[0,0],'-k')
axis('tight')
title('Business cycles (log(Y_t) - trend): HP1600 (magenta), BK(6,32,12) (green), actual cycle (blue)')

fig10 = figure;
hold
plot(timeaxis,BPLogCapitalStock,'-g')
plot(timeaxis,HPLogCapitalStock,'-m')
plot(timeaxis,LogCapitalStock-(ExpLogCapitalIntensity+LogTechnState),'-b')
plot([timeaxis(1),timeaxis(length(timeaxis))],[0,0],'-k')
axis('tight')
title('Business cycles (log(K_t)- trend): HP1600 (magenta), BK(6,32,12) (green), actual cycle (blue)')

if savepicts==1,
  if model==1,
    figure(fig1); print -depsc2 figure1a.eps  %zAggr
    figure(fig2); print -depsc2 figure1b.eps  %CapitalIntensity
    figure(fig3); print -depsc2 figure1c.eps  %OutputIntensity
    figure(fig4); print -depsc2 figure1d.eps  %LogCapitalStock
    figure(fig5); print -depsc2 figure1e.eps  %1+TechnProg
    figure(fig6); print -depsc2 figure1f.eps  %LogTechnState
    figure(fig7); print -depsc2 figure1g.eps  %Trend of log(Y_t)
    figure(fig8); print -depsc2 figure1h.eps  %Trend of log(K_t)
    figure(fig9); print -depsc2 figure1i.eps  %Business cycles (log(Y_t) - trend)
    figure(fig10); print -depsc2 figure1j.eps %Business cycles (log(K_t)- trend)
  end
  if model==2,
    figure(fig1); print -depsc2 figure2a.eps
    figure(fig2); print -depsc2 figure2b.eps
    figure(fig3); print -depsc2 figure2c.eps
    figure(fig4); print -depsc2 figure2d.eps
    figure(fig5); print -depsc2 figure2e.eps
    figure(fig6); print -depsc2 figure2f.eps
    figure(fig7); print -depsc2 figure2g.eps
    figure(fig8); print -depsc2 figure2h.eps
    figure(fig9); print -depsc2 figure2i.eps
    figure(fig10); print -depsc2 figure2j.eps
  end
  if model==3,
    figure(fig1); print -depsc2 figure3a.eps
    figure(fig2); print -depsc2 figure3b.eps
    figure(fig3); print -depsc2 figure3c.eps
    figure(fig4); print -depsc2 figure3d.eps
    figure(fig5); print -depsc2 figure3e.eps
    figure(fig6); print -depsc2 figure3f.eps
    figure(fig7); print -depsc2 figure3g.eps
    figure(fig8); print -depsc2 figure3h.eps
    figure(fig9); print -depsc2 figure3i.eps
    figure(fig10); print -depsc2 figure3j.eps
  end
end