{$R+,S+}
{$M 30000,0,200000}

(* IJIACCU.PAS

Generate the firm's time series for Salmi and Virtanen "Measuring
the Long-Run Profitability of the Firm; A Simulation Evaluation of
the Financial Statement Based IRR Estimation Methods". The version
with the double-declining balance depreciation and the estimation
with Ijiri's method with perfect knowledge of the gross book value.

Usage of the program: IJIACCU SMUCPARA.DTA

SMUCPARA.DTA contents outline:

#Identification
1980  #First year (only for computer runs; firstYear)
10    #Total length of the simulation period (T; yearsToSim)
0.08  #Growth rate (k; growthRate)
0.50  #Amplitude of the sinusoidal business cycles (A; amplitude)
0.20  #Standard deviation in the random component
2.5   #Shock coefficient (S; shock)
9999  #timing of the investment shock (ç; tau)
40    #First capital expenditure (g(0); firstCapE)
0.0   #Contribution coefficient b(0)
0.7   #Contribution coefficient b(1)
0.6   #Contribution coefficient b(2)
end   #That's it

*)

program IjiaccuProgram;

uses SMUCTPU;  { A unit of common routines }

const progDesc = 'Simulate observations; Estimate Ijiri''s accurate results';
      progName = 'IJIACCU';
      progDate = 'Fri 15-Nov-96';
      progVers = '1.0';
      maxYears = 50;

type realVectorType = array [0..maxYears] of real;
     realMatrixType = array [0..maxYears] of realVectorType;
     realMatrixPtrType = ^realMatrixType;

(* Input variables and constants *)
var firstYear,      { for illustrative purposes only }
    yearsToSim,
    tau             { timing of the investment shock 24; 30}
     : integer;

var growthRate,
    amplitude,      { amplitude of the business cycle }
    shock,          { shock coefficient }
    sigma           { standard deviation in the random component }
     : real;

const cycleLength = 6.0;  { length of the business cycle (C) }

var b               { contribution coefficients }
     : realVectorType;

(* Variables for processing *)
var p,              { accountant's profit }
    v,              { net book value }
    f,              { cash inflow }
    g,              { capital expendidure }
    d               { depreciation amount }
     : realVectorType;

                      { Accurate indicates that cumulative depreciation }
                      { is known precisely }
var accurGBV,         { accurate gross book value }
    accurD,           { accurate cumulative depreciation }
    accurCRR,         { accurate cash recovery rate }
    estimGBV,         { estimated gross book value }
    estimD,           { estimated cumulative depreciation }
    estimCRR          { estimated cash recovery rate }
     : realVectorType;

var lifeSpan,         { True life-span of the capital invesments }
    estimLifeSpan,    { Estimated life-span of the capital invesments }
    iter              { Number of iterations in caulcations }
     : integer;
var n3,               { year of last available observation }
    n2                { first CRR calculation year - 1 }
     : integer;

var true_irr,             { The true internal rate of return from }
                          { the contribution distribution }
    deprCoeff,            { double-declining depreciation coefficient }
    average_accurate_crr, { Average CRR for the observation period }
                          { with accurate culmulative depreciation known }
    accurate_irr_estim,   { Esimate of IRR from CRR with accurate cum. depr.}
    average_estim_crr,    { Average CRR with estimated accumulated depr. }
    estimated_irr_estim   { Esimate of IRR from CRR with estimated cum. depr.}
     : real;

(* Get the estimated life span from the command line *)
procedure GETSPAN (var estimLifeSpan : integer);
var s : string;
    i, k : integer;
begin
  if paramCount < 2 then begin
    writeln;
    writeln ('Usage: ', progName, ' [InputFileName EstimatedLifeSpan]');
    halt;
  end; {endif}
  s := paramStr(2);
  Val (s, estimLifeSpan, k);
  if k > 0 then begin
    writeln;
    writeln (s);
    for i := 1 to k-1 do write (' '); writeln ('^');
    writeln (#7, 'Error in integer');
    halt;
  end; {if}
end;  (* getspan *)

(* Read and display one file from the input parameter file *)
procedure READLINE (var filePointer : text;    { <-- }
                    tx              : string;  { <-- }
                    decimals        : byte;    { <-- }
                    var number      : real);   { --> }
var s : string;   { read line as string }
    k : integer;  { error index }
    i : integer;  { counter }
    j : integer;  { position of the comment initiator '#' }
begin
  i := 0;
  while (i = 0) and not eof(filePointer) do begin
    readln (filePointer, s);
    j := Pos('#',s);
    if (Length(s) > 0) and ((j = 0) or (j > 1)) then begin
      Inc (i);
      if j > 1 then s := TRAILFN (Copy (s, 1, j-1));
      Val (s, number, k);
      writeln (tx, number:0:decimals);
    end
      else writeln (s); {endif}
  end; {while}
end; (* readline *)

(* Read and display the input data *)
procedure READDATA (var firstYear  : integer;
                    var yearsToSim : integer;
                    var growthRate : real;
                    var amplitude  : real;
                    var sigma      : real;
                    var shock      : real;
                    var tau        : integer;
                    var firstCapE  : real;
                    var lifeSpan   : integer;
                    var b : realVectorType);
label endloop;
var filePointer : text;
    number      : real;       { auxiliary variable }
    s           : string;
    j, k        : integer;
begin
  {}
  {... open the input file ...}
  if ParamCount > 0 then begin
    if not FILEXIST (ParamStr(1)) then begin
      writeln;
      writeln ('File ', ParamStr(1), ' not found');
      halt;
    end; {if}
    assign (filePointer, ParamStr(1));
  end
  else halt; { this should never be reached because of getspan }
  reset (filePointer);
  {}
  {... read the parameter file contents ...}
  READLINE (filePointer, '#First year: ', 0, number);
  firstYear := Trunc(number);
  READLINE (filePointer, '#Years to simulate: ', 0, number);
  yearsToSim := Trunc(number);
  READLINE (filePointer, '#Growth rate: ', 4, growthRate);
  READLINE (filePointer, '#Amplitude of cycles: ', 4, amplitude);
  READLINE (filePointer, '#Random term STD: ', 4, sigma);
  READLINE (filePointer, '#Shock coefficient: ', 4, shock);
  READLINE (filePointer, '#Shock timing: ', 0, number);
  tau := Trunc(number);
  READLINE (filePointer, '#First capital expenditure: ',
                          4, firstCapE);
  {}
  FillChar (b, SizeOf(b), 0);
  lifeSpan := 1;
  repeat
    readln (filePointer, s);
    if eof(filePointer) then goto endloop;
    s := TRAILFN(s);
    if s = 'end' then goto endloop;
    j := Pos('#',s);
    if (Length(s) > 0) and ((j = 0) or (j > 1)) then begin
      if j > 1 then s := TRAILFN (Copy (s, 1, j-1));
      Val (s, b[lifeSpan], k);
      writeln ('#', lifeSpan : 4, b[lifeSpan]:10:4);
      Inc (lifeSpan);
    end; {if}
  until false; {repeat}
  endloop:
  lifeSpan := lifeSpan - 1;
  {}
  close (filePointer);
end;  (* readdata *)

(* Auxiliary function to calculate net present value *)
function NPVFN (b        : realVectorType;
                rate     : real;
                lifeSpan : integer) : real;
var y : real;
    i : integer;
begin
  y := -1.0;
  for i := 0 to lifeSpan do
    y := y + b[i] / GENPOWFN (1.0 + rate, i);
  npvfn := y;
end;  (* npvfn *)

(* Internal rate of return for the contribution distribution *)
function IRRFN (b : realVectorType; lifeSpan : integer) : real;
const maxIter = 30;
var x, x1, x2, y1, y2 : real;
    j : integer;
begin
  x1 := 0.0;
  x2 := 0.2;
  for j := 1 to maxIter do begin
    y1 := NPVFN (b, x1, lifeSpan);
    y2 := NPVFN (b, x2, lifeSpan);
    x := (x1*y2 - x2*y1) / (y2 - y1);
    x1 := x2;
    x2 := x;
    irrfn := x;
    if abs(x2-x1) < 0.000001 then exit;
  end; {for}
  writeln ('The IRRFN algorithm failed in ', maxIter, ' iterations');
  halt;
end;  (* irrfn *)

(* Perform the simulation *)
procedure DOSIMUL (yearsToSim  : integer;           { <-- }
                   growthRate  : real;              { <-- }
                   lifeSpan    : integer;           { <-- }
                   amplitude   : real;              { <-- }
                   sigma       : real;              { <-- }
                   shock       : real;              { <-- }
                   tau         : integer;           { <-- }
                   var g       : realVectorType;    { <--> }
                   var f,d,p,v : realVectorType);   { --> }
var t, i, m : integer;
    rn : real;
begin
  deprCoeff := 2.0 * (1.0 / lifeSpan);
  {... perform the simulation for all the years ...}
  for t := 0 to yearsToSim do begin
    if t > 0 then begin  { skip the "zero" year ...}
      rn := RNDNRMFN;
      g[t] := g[0] * GENPOWFN ((1.0 + growthRate),(t))
            * (1.0 + amplitude *
               sin(2.0*pi/cycleLength*(t)+pi/cycleLength) )
            * (1.0 + sigma*rn);
      if t = tau then begin
        g[t] := g[t] + shock * g[t];
        g[t] := g[t] /(1.0 + sigma*rn);
      end;
    end;
    {}
    m := MINIFN (lifeSpan, t);
    f[t] := 0.0;
    for i := 0 to m do
      f[t] := f[t] + b[i]*g[t-i];
    {}
    for i := 1 to m do begin
      d[t] := d[t] + deprCoeff * GENPOWFN(1.0-deprCoeff,i-1) * g[t-i];
      if (i = m) and (m = lifeSpan) then
        d[t] := d[t] + GENPOWFN(1.0-deprCoeff,m) * g[t-m];
    end; {for i}
    {}
    p[t] := f[t] - d[t];
    {}
    if t > 0 then v[t] := v[t-1] + g[t] - d[t]
      else v[t] := g[t];
  end; {for t}
end;  (* dosimul *)

(* Calculate the depreciation matrix, accurate book value, annual CRRs,
   with declining balance depreciation *)
procedure CALC_ACCUR_GROSS_BOOK_VALUE_ETC
                                 (yearsToSim   : integer;
                                  lifeSpan     : integer;
                                  g            : realVectorType;
                                  var accurGBV : realVectorType;
                                  var accurD   : realVectorType;
                                  var accurCRR : realVectorType);
var dmPtr   { declining depreciation for each capital investment and year }
      : realMatrixPtrType;
var t, i, j : integer;
    sum : real;
    deprCoeff : real;   { Double declining balance depreciation coefficient }
begin
  {... Calculate double declining balance depreciation ...}
  deprCoeff := 2.0 * (1.0 / lifeSpan);
  New (dmPtr);
  FillChar (dmPtr^, SizeOf(realMatrixType), 0);
  for i := 0 to yearsToSim-1 do begin
    for j := i+1 to MINIFN (i+lifeSpan, yearsToSim) do begin
      { declining balance depreciation }
      dmPtr^[i,j] := deprCoeff * GENPOWFN(1.0-deprCoeff,j-i-1) * g[i];
      if (j = i+lifeSpan) then
        dmPtr^[i,j] := dmPtr^[i,j]
                     + GENPOWFN(1.0-deprCoeff,lifeSpan) * g[i];
    end; {for j}
  end; {for i}
  {}
  {... Gross capital investments V(t) = v(t)+cumulative depr ...}
  {... initialize vectors ...}
  FillChar (accurGBV, SizeOf(realVectorType), 0);
  FillChar (accurD, SizeOf(realVectorType), 0);
  FillChar (accurCRR, SizeOf(realVectorType), 0);
  accurGBV[0] := v[0];
  for t := 1 to yearsToSim do begin
    sum := 0.0;
    for i := MAXIFN (0, t+1-lifeSpan) to t-1 do
      for j := MAXIFN (1, t-lifeSpan+1) to t do
        sum := sum + dmPtr^[i,j];
    accurGBV[t] := v[t] + sum;
    accurD[t] := sum;
    accurCRR[t] := f[t] / accurGBV[t-1];
  end; {for t}
  {}
  Dispose (dmPtr);
end;  (* calc_accur_gross_book_value_etc *)

(* Calculate the depreciation matrix, estimated groos book value,
   annual CRRs *)
procedure CALC_ESTIM_GROSS_BOOK_VALUE_ETC
                                 (yearsToSim    : integer;
                                  trueLifeSpan  : integer;
                                  estimLifeSpan : integer;
                                  f             : realVectorType;
                                  d             : realVectorType;
                                  v             : realVectorType;
                                  var estimGBV  : realVectorType;
                                  var estimD    : realVectorType;
                                  var estimCRR  : realVectorType;
                                  var n2, n3    : integer);
var t, i : integer;
    sum : real;
    accuYears : integer;   { depreciation accumulation span }
    n1 : integer;          { year of first available observation }
begin
  {}
  {... Gross capital investments V(t) = v(t)+cumulative depr ...}
  {... initialize vectors ...}
  FillChar (estimGBV, SizeOf(realVectorType), 0);
  FillChar (estimD, SizeOf(realVectorType), 0);
  FillChar (estimCRR, SizeOf(realVectorType), 0);
  {}
  accuYears := (estimLifeSpan+1) div 2;   {The +1 is to tackle odd life-spans}
  n3 := yearsToSim;           { year of last available observation }
  n1 := trueLifeSpan + 1;     { year of first available observation }
  n2 := n1 + accuYears;       { first CRR calculation year - 1 }
  for t := n3 downto n2 do begin
    estimD[t] := 0; { The estimate of the accumulated depreciation }
    for i := t downto t-accuYears+1 do begin
      estimD[t] := estimD[t] + d[i];
    end; {for i}
    estimGBV[t] := v[t] + estimD[t];
  end; {for t}
  {}
  for t := n3 downto n2+1 do begin
    estimCRR[t] := f[t] / estimGBV[t-1];
  end; {for t}
  {}
  sum := 0;
  for t := n3 downto n2+1 do begin
    sum := sum + estimCRR[t];
  end; {for t}
end;  (* calc_estim_gross_book_value_etc *)

function AVERAGE_CRR_FN (crr : realVectorType;
                         m1  : integer;
                         m2  : integer) : real;
var t : integer;
    sum : real;
begin
  sum := 0;
  for t := m1 to m2 do sum := sum + crr[t];
  average_crr_fn := sum / (m2-m1+1);
end;  (* average_crr_fn *)

(* Ijiri's estimation formula *)
function FFN (irr, crr : real; n : integer) : real;
begin
  ffn := crr - irr/(1.0-GENPOWFN((1.0+irr),-n));
end;  (* fnn *)

function IRRESTFN (crr : real;
                   n : integer;
                   var iter : integer) : real;
const maxIter = 30;
var a, a1, a2, f1, f2 : real;
    i : integer;
begin
  a1 := 0.1;
  a2 := 0.2;
  for i := 1 to maxIter do begin
    f1 := FFN (a1, crr, n);
    f2 := FFN (a2, crr, n);
    a  := (a1 * f2 - a2 * f1) / (f2 - f1);
    a1 := a2;
    a2 := a;
    irrestfn := a2;
    iter := i;
    if (abs(a2-a1) < 0.00001) then exit;
  end; {for}
  writeln ('The IRRESTFN algorithm failed in ', maxIter, ' iterations');
  halt;
end;  (* irrestfn *)

(* Write the simulated data *)
procedure WRITEDATA (firstYear  : integer;
                     yearsToSim : integer;
                     lifeSpan   : integer;
                     g,f,d,p,v  : realVectorType;
                     irr        : real);
const ff1 =  5;
      ff2 = 12;
      dd  =  4;
var t : integer;
begin
  writeln;
  write ('#Year'      : ff1);
  write ('Capital'    : ff2);
  write ('FundsFrom'  : ff2);
  write ('Ddeclining' : ff2);
  write ('Operating'  : ff2);
  write ('Net book'   : ff2);
  writeln;
  {}
  write ('#', ' '     : ff1-1);
  write ('expendit'   : ff2);
  write ('operations' : ff2);
  write ('depreciat'  : ff2);
  write ('income'     : ff2);
  write ('value'      : ff2);
  writeln;
  {}
  write ('#', ' '     : ff1-1);
  write ('g(t)'       : ff2);
  write ('f(t)'       : ff2);
  write ('d(t)'       : ff2);
  write ('p(t)'       : ff2);
  write ('v(t)'       : ff2);
  writeln;
  {}
  for t := 0 to yearsToSim do begin
    if t <= lifeSpan + 1 then write ('#', firstYear + t : ff1-1)
      else write (firstYear + t : ff1);
    write (g[t] : ff2:dd);
    write (f[t] : ff2:dd);
    write (d[t] : ff2:dd);
    write (p[t] : ff2:dd);
    write (v[t] : ff2:dd);
    writeln;
  end; {for}
  {}
  writeln;
  writeln ('#True internal rate of return = ', 100.0*irr:0:4, '%');
end;  (* writedata *)

(* Write the simulated data *)
procedure WRITEDATA2 (firstYear            : integer;
                      yearsToSim           : integer;
                      lifeSpan             : integer;
                      accurGBV             : realVectorType;
                      accurD               : realVectorType;
                      accurCRR             : realVectorType;
                      estimGBV             : realVectorType;
                      estimD               : realVectorType;
                      estimCRR             : realVectorType;
                      average_accurate_crr : real;
                      accurate_irr_estim   : real;
                      average_estim_crr    : real;
                      estimated_irr_estim  : real;
                      irr                  : real);
const ff1 =  5;
      ff2 = 12;
      dd  =  4;
var t : integer;
begin
  writeln;
  write ('#Year'      : ff1);
  write ('AccurCumul' : ff2);
  write ('AccurGross' : ff2);
  write ('Accur CRR'  : ff2);
  write ('Est. cumul' : ff2);
  write ('Est. gross' : ff2);
  write ('Est. CRR'   : ff2);
  writeln;
  {}
  write ('#', ' '     : ff1-1);
  write ('depreciat'  : ff2);
  write ('book value' : ff2);
  write (''           : ff2);
  write ('depreciat'  : ff2);
  write ('book value' : ff2);
  write (''           : ff2);
  writeln;
  {}
  write ('#', ' '     : ff1-1);
  write ('AccurD(t)'  : ff2);
  write ('AccurV(t)'  : ff2);
  write ('AccuCRR(t)' : ff2);
  write ('D(t)'       : ff2);
  write ('V(t)'       : ff2);
  write ('CRR(t)'     : ff2);
  writeln;
  {}
  for t := 0 to yearsToSim do begin
    if t <= lifeSpan + 1 then write ('#', firstYear + t : ff1-1)
      else write (firstYear + t : ff1);
    write (accurD[t]   : ff2:dd);
    write (accurGBV[t] : ff2:dd);
    write (accurCRR[t] : ff2:6);
    write (estimD[t]   : ff2:dd);
    write (estimGBV[t] : ff2:dd);
    write (estimCRR[t] : ff2:6);
    writeln;
  end; {for}
  {}
  writeln;
  writeln ('#Average accurate CRR           = ',
            100*average_accurate_crr:ff2:dd, '%');
  writeln ('#IRR from average accurate CRR  = ',
            100*accurate_irr_estim:ff2:dd, '%');
  writeln;
  writeln ('#Average estimated CRR          = ',
            100*average_estim_crr:ff2:dd, '%');
  writeln ('#IRR from average estimated CRR = ',
            100*estimated_irr_estim:ff2:dd, '%');
  writeln;
  writeln ('#True internal rate of return   = ',
            100.0*irr:ff2:dd, '%');
end;  (* writedata2 *)

(* Main program *)
begin
  FileMode := 0;
  HEADER (progName, progDesc, progVers, progDate, '#');
  GETSPAN (estimLifeSpan);
  READDATA (firstYear, yearsToSim, growthRate,
            amplitude, sigma, shock, tau, g[0], lifeSpan, b);
  true_irr := IRRFN (b, lifeSpan);
  DOSIMUL (yearsToSim, growthRate, lifeSpan, amplitude, sigma,
           shock, tau, g, f, d, p, v);
  WRITEDATA (firstYear, yearsToSim, lifeSpan, g, f, d, p, v, true_irr);
  writeln;
  {}
  CALC_ESTIM_GROSS_BOOK_VALUE_ETC
    (yearsToSim, lifeSpan, estimLifeSpan, f, d, v,    { <-- }
     estimGBV, estimD, estimCRR, n2, n3);             { --> }
  average_estim_crr := AVERAGE_CRR_FN (estimCRR, n2+1, n3);
  estimated_irr_estim := IRRESTFN (average_estim_crr, estimLifeSpan, iter);
  {}
  CALC_ACCUR_GROSS_BOOK_VALUE_ETC (yearsToSim, lifeSpan, g,
                                   accurGBV, accurD, accurCRR);
  average_accurate_crr := AVERAGE_CRR_FN (accurCRR, n2+1, n3);
  accurate_irr_estim := IRRESTFN (average_accurate_crr,
                                  estimLifeSpan, iter);
  {}
  WRITEDATA2 (firstYear, yearsToSim, lifeSpan,
              accurGBV, accurD, accurCRR,
              estimGBV, estimD, estimCRR,
              average_accurate_crr,
              accurate_irr_estim,
              average_estim_crr,
              estimated_irr_estim,
              true_irr);
end.  (* IjiaccuProgram *)