function run_model(simcase);
% function run_model(case);  case = 1 or 2
%    Simulates the OLG model introduced in the paper
%    Economic Growth and Business Cycles: A Critical Comment on Detrending Time Series
%    by Klaus Reiner Schenk-Hoppé
%    Working Paper, August 2000
%    Institute for Empirical Research in Economics, University of Zurich, Switzerland
%
%    case 1: parameters defined in Section 4.1
%    case 2: parameters defined in Section 4.2
%
% Author: Klaus Reiner Schenk-Hoppé (klaus@iew.unizh.ch)
% Last update: August 15, 2000


global alpha_RBC
global beta_RBC
global gamma_RBC
global delta_RBC
global L_RBC
global w_RBC
global ExpLaz_RBC

alpha_RBC = 0.42;
beta_RBC = 0.9883;
delta_RBC = 0.02599;
gamma_RBC = 3.0;

M=208;
timeaxis = (0:0.25:M*4-0.25);

ORIGL_RBC(1) = 1.0;

for t = 1:M,
  ORIGL_RBC(t+1) = 1.00 * ORIGL_RBC(t);
end;

L_RBC = ORIGL_RBC(1);

A = 0.95;
rho2 = 0.005;
rho1 = - rho2;
xi_s = 0.0;
epsilon = 0.000001;
z_RBC = 1.0;
Ez_RBC = 1.0;

if (simcase == 1),
  B = 0.95;
  phi2 = 0.005;
else
  B = 0.99;
  phi2 = 0.0005;
end;

C = 0.02;
phi1 = - phi2;
eta_s = 0.02;

a_RBC = 1.0;
k(1) = 10000000;
a(1) = a_RBC;
z(1) = Ez_RBC;
xi(1) = xi_s;
eta(1) = eta_s;

%check validity of parameters
if (phi2 >= (1-B)*(1+C)),
  fprintf('Condition (ii) of Assumption 4.1 is violated \n')
  return
end

rst = rand('state');

%%%%%%%%%%%% pre-iteration
for t = 1:50,
  L_RBC = ORIGL_RBC(1);
  z_RBC = z(1);
  a_RBC = a(1);

  wage(t) = (1-alpha_RBC) .* z(1) .* a(1) .* (k(t)./(a(1).*ORIGL_RBC(1))).^alpha_RBC;
  interest(t) = alpha_RBC .* z(1) .* ((a(1).*ORIGL_RBC(1))./k(t)).^(1-alpha_RBC);
  output(t) = z(1) .* k(t).^alpha_RBC .* (a(1).*ORIGL_RBC(1)).^(1-alpha_RBC);

  ExpLaz_RBC = (ORIGL_RBC(1) * a(1)).^(1.0-alpha_RBC) * 1.0/(2.0-alpha_RBC) *((1+B*eta(1)+phi2).^(2.0-alpha_RBC) - (1+B*eta(1)+phi1).^(2.0-alpha_RBC))*( 1.0 + 2.0*(1.0-epsilon)/pi *( (A*xi(1)+rho2)*atan(A*xi(1)+rho2) - 0.5*log(1.0+(A*xi(1)+rho2)^2.0)-(A*xi(1)+rho1)*atan(A*xi(1)+rho1) + 0.5*log(1.0+(A*xi(1)+rho1)^2.0)));

  w_RBC = wage(t);
  mu(t) = fzero('ree',3);

  s(t) = wage(t)/(1.0 + (beta_RBC* (1-delta_RBC + mu(t))^(-1.0/gamma_RBC)));

  k(t+1) = ORIGL_RBC(1) * s(t);
end;

k(1) = k(50);
a(1) = a_RBC;
z(1) = Ez_RBC;
xi(1) = xi_s;
eta(1) = eta_s;

for t = 1:M,
  L_RBC = ORIGL_RBC(t);
  z_RBC = z(t);
  a_RBC = a(t);

  wage(t) = (1-alpha_RBC) .* z(t) .* a(t) .* (k(t)./(a(t).*ORIGL_RBC(t))).^alpha_RBC;
  interest(t) = alpha_RBC .* z(t) .* ((a(t).*ORIGL_RBC(t))./k(t)).^(1-alpha_RBC);
  output(t) = z(t) .* k(t).^alpha_RBC .* (a(t).*ORIGL_RBC(t)).^(1-alpha_RBC);

  ExpLaz_RBC = (ORIGL_RBC(t+1) * a(t)).^(1.0-alpha_RBC) * 1.0/(2.0-alpha_RBC) *((1+B*eta(t)+phi2).^(2.0-alpha_RBC) - (1+B*eta(t)+phi1).^(2.0-alpha_RBC))*( 1.0 + 2.0*(1.0-epsilon)/pi *( (A*xi(t)+rho2)*atan(A*xi(t)+rho2) - 0.5*log(1.0+(A*xi(t)+rho2)^2.0)-(A*xi(t)+rho1)*atan(A*xi(t)+rho1) + 0.5*log(1.0+(A*xi(t)+rho1)^2.0)));

  w_RBC = wage(t);
  mu(t) = fzero('ree',3);

  s(t) = wage(t)/(1.0 + (beta_RBC* (1-delta_RBC + mu(t))^(-1.0/gamma_RBC)));

  k(t+1) = ORIGL_RBC(t) * s(t);

  eta(t+1) = B * eta(t) + (1-B)*C + phi1 + (phi2-phi1)*rand;
  a(t+1) = (1+ eta(t+1)) * a(t);

  xi(t+1) = A * xi(t) + rho1 + (rho2-rho1)*rand;
  z(t+1) = 1+2.0*(1-epsilon)/pi * atan(xi(t+1));
end;

%%%% quarterly data adjusted at annual data
TotalOutput = 4.0*output(1:M);
TotalOutputAdj = 4.0*output(1:M)./z(1:M); % adjusted for total hours worked
NetInterestRate = interest(1:M).^4;
GrossInterestRate = (1-delta_RBC+interest(1:M)).^4;
CapitalStock = 4.0.*k(1:M);
CapitalPerWorker = 4.0.*k(1:M)./ORIGL_RBC(1:M);
WageIncome = 4.0.*wage(1:M);
WageRate = 4.0.*wage(1:M)./z(1:M);  % wage rate per unit of labor (i.e. adjusted for total hours worked)
OutputPerWorker = TotalOutput(1:M)./ORIGL_RBC(1:M);
Output = TotalOutput(1:M)./(z(1:M).*ORIGL_RBC(1:M)); % per unit of labor (i.e. adjusted for total hours worked)
TotalWageIncome = WageIncome .* ORIGL_RBC(1:M);
TotalCapitalIncome = (1-delta_RBC).^4 .* CapitalStock + TotalOutput - TotalWageIncome;


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

sGDP = hpfilter(log(TotalOutput),1600);
HPGDP = sGDP*transpose(log(TotalOutput));
HPGDPAdj = hpfilter(log(TotalOutputAdj),1600)*transpose(log(TotalOutputAdj));
HPCapitalStock = hpfilter(log(CapitalStock),1600)*transpose(log(CapitalStock));
HPCapitalPerWorker = hpfilter(log(CapitalPerWorker),1600)*transpose(log(CapitalPerWorker));
HPWageIncome = hpfilter(log(WageIncome),1600)*transpose(log(WageIncome));
HPWageRate = hpfilter(log(WageRate),1600)*transpose(log(WageRate));
HPOutputPerWorker = hpfilter(log(OutputPerWorker),1600)*transpose(log(OutputPerWorker));
HPOutput = hpfilter(log(Output),1600)*transpose(log(Output));
HPTechnProg = hpfilter(log(a),1600)*transpose(log(a));

HPtrendGDP = log(TotalOutput(1:M))-transpose(HPGDP(1:M));
HPtrendGDPAdj = log(TotalOutputAdj(1:M))-transpose(HPGDPAdj(1:M));
HPtrendCapitalStock = log(CapitalStock(1:M))-transpose(HPCapitalStock(1:M));
HPtrendCapitalPerWorker = log(CapitalPerWorker(1:M))-transpose(HPCapitalPerWorker(1:M));
HPtrendWageIncome = log(WageIncome(1:M))-transpose(HPWageIncome(1:M));
HPtrendWageRate = log(WageRate(1:M))-transpose(HPWageRate(1:M));
HPtrendOutputPerWorker = log(OutputPerWorker(1:M))-transpose(HPOutputPerWorker(1:M));
HPtrendOutput = log(Output(1:M))-transpose(HPOutput(1:M));
HPtrendTechnProg = log(a(1:M))-transpose(HPTechnProg(1:M));

HPWageIncome4 = hpfilter(log(WageIncome),4)*transpose(log(WageIncome));
HPWageIncome40 = hpfilter(log(WageIncome),40)*transpose(log(WageIncome));
HPWageIncome400 = hpfilter(log(WageIncome),400)*transpose(log(WageIncome));

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

BPGDP = bpf(transpose(log(TotalOutput)),6,32,12);
BPGDPAdj = bpf(transpose(log(TotalOutputAdj)),6,32,12);
BPCapitalStock = bpf(transpose(log(CapitalStock)),6,32,12);
BPCapitalPerWorker = bpf(transpose(log(CapitalPerWorker)),6,32,12);
BPWageIncome = bpf(transpose(log(WageIncome)),6,32,12);
BPWageRate = bpf(transpose(log(WageRate)),6,32,12);
BPOutputPerWorker = bpf(transpose(log(OutputPerWorker)),6,32,12);
BPOutput = bpf(transpose(log(Output)),6,32,12);
BPTechnProg = bpf(transpose(log(a)),6,32,12);

BPtrendGDP = log(TotalOutput(1:M))-transpose(BPGDP(1:M));
BPtrendGDPAdj = log(TotalOutputAdj(1:M))-transpose(BPGDPAdj(1:M));
BPtrendCapitalStock = log(CapitalStock(1:M))-transpose(BPCapitalStock(1:M));
BPtrendCapitalPerWorker = log(CapitalPerWorker(1:M))-transpose(BPCapitalPerWorker(1:M));
BPtrendWageIncome = log(WageIncome(1:M))-transpose(BPWageIncome(1:M));
BPtrendWageRate = log(WageRate(1:M))-transpose(BPWageRate(1:M));
BPtrendOutputPerWorker = log(OutputPerWorker(1:M))-transpose(BPOutputPerWorker(1:M));
BPtrendOutput = log(Output(1:M))-transpose(BPOutput(1:M));
BPtrendTechnProg = log(a(1:M))-transpose(BPTechnProg(1:M));



%%%%%%%%%%%%%% linear trend

[lr_GDP,lr_dataLogGDP] = lin_regress(timeaxis(1:M),log(TotalOutput(1:M)));
[lr_GDPAdj,lr_dataLogGDPAdj] = lin_regress(timeaxis(1:M),log(TotalOutputAdj(1:M)));
[lr_CapitalStock,lr_dataLogCapitalStock] = lin_regress(timeaxis(1:M),log(CapitalStock(1:M)));
[lr_CapitalPerWorker,lr_dataLogCapitalPerWorker] = lin_regress(timeaxis(1:M),log(CapitalPerWorker(1:M)));
[lr_WageIncome,lr_dataLogWageIncome] = lin_regress(timeaxis(1:M),log(WageIncome(1:M)));
[lr_WageRate,lr_dataLogWageRate] = lin_regress(timeaxis(1:M),log(WageRate(1:M)));
[lr_OutputPerWorker,lr_dataLogOutputPerWorker] = lin_regress(timeaxis(1:M),log(OutputPerWorker(1:M)));
[lr_Output,lr_dataLogOutput] = lin_regress(timeaxis(1:M),log(Output(1:M)));
[lr_TechnProg,lr_dataLogTechnProg] = lin_regress(timeaxis(1:M),log(a(1:M)));

l_detrended_LogGDP = log(TotalOutput) - lr_dataLogGDP;
l_detrended_LogGDPAdj = log(TotalOutputAdj) - lr_dataLogGDPAdj;
l_detrended_LogCapitalStock = log(CapitalStock) - lr_dataLogCapitalStock;
l_detrended_LogCapitalPerWorker = log(CapitalPerWorker) - lr_dataLogCapitalPerWorker;
l_detrended_LogWageIncome = log(WageIncome) - lr_dataLogWageIncome;
l_detrended_LogWageRate = log(WageRate) - lr_dataLogWageRate;
l_detrended_LogOutputPerWorker = log(OutputPerWorker) - lr_dataLogOutputPerWorker;
l_detrended_LogOutput = log(Output) - lr_dataLogOutput;
l_detrended_LogTechnProg = log(a(1:M)) - lr_dataLogTechnProg;


if (simcase == 1),
   save case1data;
else
   save case2data;
end;

fig1 = figure;
hold
plot(timeaxis(1:M),z(1:M),'-b');
plot([timeaxis(1),timeaxis(M)],[1,1],'-k')
grid on
axis('tight')
title('Total units of labor employed z_t. The time scale is in years and divided into quarters.');
%print -depsc2 figure1.eps

fig2 = figure;
hold
plot(timeaxis(1:M),1+eta(1:M),'-b');
plot([timeaxis(1),timeaxis(M)],[1.02,1.02],'-k')
grid on
axis('tight')
title('Rate of labor-augmenting technical progress 1+eta_{t+1}.');
%print -depsc2 figure2.eps

fig3 = figure;
hold
plot(timeaxis(1:M),l_detrended_LogGDP(1:M),'-b')
plot(timeaxis(1:M),HPtrendGDP(1:M) - lr_dataLogGDP(1:M),'-g')
plot(timeaxis(1:M),BPtrendGDP(1:M) - lr_dataLogGDP(1:M),'-m')
plot(timeaxis(1:M),l_detrended_LogTechnProg,'-r')
plot([timeaxis(1),timeaxis(M)],[0,0],'-k')
grid on
axis('tight')
title('Deviation from linear trend of log(Y_t): data (blue), HP(1600) trend (green), BP(6,32,12) trend (magenta), and linear detrended log(a_t) (red)')
%print -depsc2 figure3.eps


fig4 = figure;
hold
plot(timeaxis(1:M),HPGDP(1:M),'-g')
plot(timeaxis(1:M),BPGDP(1:M),'-m')
plot(timeaxis(1:M),l_detrended_LogGDP(1:M)-l_detrended_LogTechnProg,'-b')
plot([timeaxis(1),timeaxis(M)],[0,0],'-k')
axis('tight')
title('Business cycles of GDP (deviation of log(Y_t) from different trends): HP(1600) trend (green), BP(6,32,12) trend (magenta), sum of the linear detrended log(a_t) data and linear trend of log(Y_t) (blue)')
%print -depsc2 figure4.eps

time1 = 27*4;
time2 = M;
%time1 = 1;
%time2 = M;
fig5 = figure;
hold
plot(timeaxis(time1:time2),HPWageIncome4(time1:time2),'-r')
plot(timeaxis(time1:time2),HPWageIncome40(time1:time2),'-g')
plot(timeaxis(time1:time2),HPWageIncome400(time1:time2),'-c')
plot(timeaxis(time1:time2),BPWageIncome(time1:time2),'-m')
plot(timeaxis(time1:time2),l_detrended_LogWageIncome(time1:time2)-l_detrended_LogTechnProg(time1:time2),'-b')
plot([timeaxis(time1),timeaxis(time2)],[0,0],'-k')
axis('tight')
title('Business cycles of wage income (deviation of log(w_t) from different trends): HP(4) trend (red), HP(40) trend (green), HP(400) trend (cyan), BP(6,32,12) trend, sum of the linear detrended log(a_t) data and linear trend of log(w_t)')
%print -depsc2 figure5.eps

HP200 = hpfilter(log(TotalOutput(1:M)),200)*transpose(log(TotalOutput(1:M)));
HP3200 = hpfilter(log(TotalOutput(1:M)),3200)*transpose(log(TotalOutput(1:M)));
HP6400 = hpfilter(log(TotalOutput(1:M)),6400)*transpose(log(TotalOutput(1:M)));

fig5c = figure;
hold
plot(timeaxis(1:M),l_detrended_LogGDP(1:M)-l_detrended_LogTechnProg,'-b');
plot(timeaxis(1:M),HP6400(1:M),'-m');
plot(timeaxis(1:M),HP3200(1:M),'-g');
plot(timeaxis(1:M),HP200(1:M),'-r');
grid on
axis('tight')
title('Business cycles of GDP (deviation of log(Y_t) from different trends): HP(6400) trend (magenta), HP(3200) trend (green), HP(200) trend (red), sum of the linear detrended log(a_t) data and linear trend of log(Y_t) (blue).');
%print -depsc2 figure5c.eps