C  SIMULATION - VALUATION OF CAPACITY INVESTMENT PROBLEM
C  WHEN SUBJECT TO A PRICE CAP
C  IN WHICH THERE IS AN INITIAL CHOICE OF CAPACITY FOLLOWED BY
C  FURTHER ADDITIONS TO CAPACITY.
C  This programme uses NAG subroutine library
C  Compiles under Fortran compilers such as f77 or f90
C  Working program but needs more development/annotation
C  Uses "in" as input file to read parameter values
C  Produces "Out" as output file
C  Outputs the expected net present value of the firm
C  based on the input parameter values read from "in"
C  Value is computed as the average value over a number of simulation paths
C  The consequence of setting the price cap at, or above, or below
C  its optimal value is considered, and for each of these scenarios
C  the value consequence of the firm setting a trigger price at, or above, 
C  or below its optimal level is then assessed
C  A total of 9 scenarios.  
C  AK=unit capacity cost at time 0 
C  AKT is the value of this at time t 
C  G=demand elasticity
C  ETA=inverse demand elasticity = 1/G 
C  ETA2 is ETA squared, VAR is SIG squared 
C  AT=level of demand at time t 
C  A0 is base level of demand in demand function at time 0
C  QT=AT*(PT^G) 
C  QT, PT are capacity, price at time t
C  PV is the trigger price
C  UP is the average rate of increase of price
C  ANPV is net present value
C  PHAT is test price (test if above or below PV)
C  EPHAT is relative test price
C  AINVT is level of investment at time t 
C  AL is alpha, the rate of growth of demand 
C  R is discount rate (annual, continuous compounded) 
C  THETA is physical depreciation rate (ditto) 
C  SIG - volatility (s.d.) of demand 
C  R - rate of interest continuous time annual
C  EPS=random error N(0,1)
C  AINVT - investment at time t
C  EPC,EPU,EPM,EPCAP are relative trigger prices 
C  under certainty, uncertainty, monopoly and price cap 
C  monopoly. EPBAR is the price cap (relative) 
C  PC,PU,PM,PCAP are associated absolute trigger prices 
C  PBAR is price cap
C  DEL is the rate of technical progress
C  ALAN1,ALAN2 landa1 landa2, the characteristic roots
C  AMUL = option multiplier
C  XP(i) is used as an arbitrary relative price at which 
C  investment is triggered.  It is set at the optimal  
C  value and at values above and below - to check that 
C  value is greatest when the trigger price is set at EPCAP
C  EPBAR is also varied above and below the EPU value
C********************************************
C  Operating variables and routines
C  NT - number of years over which value
C       is to be calculated
C  DELTA - number of subperiods
C  NTD - total number of periods =NT/DELTA
C  NREPS - number of repetitions of the simulation path
C  NPARS - number of parameter sets to be considered
C         (may not be necessary)
C   G05DDF(a,b) normal with mean a, sig b
C 
C*******************************************  
C******************************************** 
      INTEGER NT,NTD,NREPS  
      INTEGER I1,I2,I3,I4              

      DOUBLE PRECISION XP(20)                                             
      DOUBLE PRECISION AA1,AA2,AA3,AA4,DEL,AKT
      DOUBLE PRECISION EPHAT
      DOUBLE PRECISION EPU,EPC,EPM,EPBAR,EPCAP
      DOUBLE PRECISION DELTA,AL,PT,PTHAT    
      DOUBLE PRECISION G05DDF,G,AK,ETA2,VAR,UP,R1,R2
      DOUBLE PRECISION ALAN1,ALAN2,AMUL,DIFF0 
      DOUBLE PRECISION ALOGT,SIGDEL,PROFT,PVPROF,AINVT 
      DOUBLE PRECISION R,THETA,SIG,A0,ETA,DIFF1,DIFF2
      DOUBLE PRECISION AT,QT,QTHAT,EPS,VNPV,SDENPV
      DOUBLE PRECISION SDNPV,ANPV,SNPV,SSNPV,ENPV 
      EXTERNAL G05DDF,G05CBF 
C**************************************************** 
C     INPUT  section  
C**************************************************** 
      OPEN(UNIT=1,FILE='out')  
      OPEN(UNIT=2,FILE='in')  
      READ(2,8)NT,NREPS  
    8 FORMAT(3I10)  
      WRITE(6,9)NT,NREPS  
      WRITE(1,9)NT,NREPS  
    9 FORMAT('NT=Simulation path length, in years=     ',I8
     1/'NREPS=Number of paths in calculating ENPV=',I8)
C**************************************************** 
C* initialisation  
C*********************************
      CALL G05CBF(0) 
C*************************************************** 

C***************************************************  
C SET INITIAL PARAMETER VALUES 
C************************************************** 
      READ(2,11)R,THETA,AL,G,SIG,A0,AK,DEL,DELTA
   11 FORMAT(9F10.4) 
      WRITE(1,13)R,THETA,AL,G,SIG,A0,AK,DEL,DELTA
   13 FORMAT('Parameters'/
     1'r=    ',F8.4/
     1'theta=',F8.4/
     1'alpha=',F8.4/
     1'gamma=',F8.4/
     1'sigma=',F8.4/
     1'A0=   ',F8.4/
     1'K=    ',F8.4/
     1'Del,t.p.=',F8.4/
     1'Delta,sub-year time unit=',F8.4)
       ETA=1.0D0/G 
C* calculation of trigger prices  
      ETA2=ETA*ETA  
      VAR=SIG*SIG  
      UP=-ETA*(AL+THETA-(0.5D0*(ETA+1.0D0)*VAR))
      R1=UP+DEL-(0.5D0*ETA2*VAR) 
      R2=(R1*R1+(2.0D0*ETA2*VAR*(THETA+R+DEL)))**0.5D0 
      ALAN1=(-R1+R2)/(ETA2*VAR)  
      ALAN2=-(R1+R2)/(ETA2*VAR)  
      AMUL=ALAN1/(ALAN1-1.0D0)  
      EPC=THETA+R+DEL 
      EPU=AMUL*(THETA+R-UP)  
      EPM=(G/(1.0D0+G))*EPU  
      WRITE(1,15)EPC,EPU,EPM  
   15 FORMAT('Certainty relative investment trigger
     1 price,EPC=',F7.3
     1/'Uncertainty relative investment trigger price,
     1EPU=     ',F7.3
     1/'Uncertainty Monopoly relative investment
     1 trigger price,EPM=',F7.3)
      AA1=EPU-EPC 
      AA2=1.0D0+(ETA*ALAN2) 
      EPBAR=EPU-DIFF0 
      DIFF0=(EPU-EPC)/2.0D0 
      DIFF1=(EPM-EPU)/2.0D0 

C**********************************  
C** do loop over the number of VARIATIONS IN pbar  
C***********************************  
      SIGDEL=SIG*(DELTA**0.5D0)
      EPBAR=EPU-DIFF0
      DO 1100 I4=0,2
      EPBAR=EPBAR+(DIFF0*I4)
      AA3=EPBAR-EPC
      AA4=EPU*(EPBAR**(ALAN2-1.0D0))
      EPCAP=((AA3*AA4)/(AA1*AA2))**(1.0D0/ALAN2)
      WRITE(1,200)EPBAR,EPCAP
  200 FORMAT('Setting EPBAR at ',F8.4,/   
     1'this means EPCAP=',F8.4)
      DIFF2=(EPCAP-EPU)/2.0D0
      XP(1)=EPCAP-DIFF1 
      XP(2)=EPCAP   
      XP(3)=EPCAP+DIFF2 
C*************************** 
C DO LOOP OVER PARAMETER SETS   
      DO 1000 I1=1,3 
C* initialise 
      SNPV=0.0D0  
      VNPV=0.0D0  
      ENPV=0.0D0 
      SDENPV=0.0D0
      SDNPV=0.0D0 
C***************************  
C   begin do loop for the number of simulation runs  
C***************************  
      DO 900 I2=1,NREPS 
      AKT=AK
      QT=A0*((XP(I1)*AKT)**G)  
      ANPV=-AK*QT 
      ALOGT=DLOG(A0)
C****************************************  
C   do loop over number of time periods  
C****************************************  
      NTD=NT/DELTA 
      DO 800 I3=1,NTD
C at time t 
C  first sample error         
      EPS=G05DDF(0.0D0,1.0D0)  
      ALOGT=ALOGT+AL*DELTA+SIGDEL*EPS  
      AT=DEXP(ALOGT) 
C  capacity if no investment would be 
      QTHAT=QT*(DEXP(-THETA*DELTA))  
C Price of new equipment is 
      AKT=AKT*(DEXP(-DEL*DELTA)) 
C the free demand price would be 
      PTHAT=(AT**(-ETA))*(QTHAT**ETA)        
      EPHAT=PTHAT/AKT 
C If PTHAT is above the trigger price then there is investment
      IF(EPHAT.GE.XP(I1))THEN  
            QT=AT*((AKT*XP(I1))**G)  
            PT=EPBAR*AKT  
            AINVT=QT-QTHAT
      ELSE IF(EPHAT.LT.EPCAP)THEN  
            QT=QTHAT
            PT=PTHAT
            AINVT=0.0D0
      ELSE
            QT=AT*((AKT*XP(I1))**G)  
            PT=EPBAR*AKT   
            AINVT=0.0D0  
      ENDIF

      CONTINUE 
      PROFT=DELTA*PT*QT-AKT*AINVT  
      PVPROF=PROFT*DEXP(-R*DELTA*I3)   
      ANPV=ANPV+PVPROF 
  800 ENDDO 
C**************************************           
      SNPV=SNPV+ANPV 
      SSNPV=SSNPV+ANPV*ANPV   
  900 ENDDO 
C**************************************            
      ENPV=SNPV/NREPS  
      VNPV=((SSNPV/NREPS)-(ENPV*ENPV))*NREPS/(NREPS-1) 
      IF(VNPV.LE.0.0D0)THEN   
         SDNPV=0.0D0  
         GOTO 820  
      ENDIF 
      SDNPV=VNPV**0.5D0  
  820 CONTINUE
      SDENPV=SDNPV/(NREPS**0.5D0)  
      WRITE(1,850)
      WRITE(1,853)XP(I1)
      WRITE(1,851)ENPV,SDENPV
  850 FORMAT(10(F12.4,1X)) 
  851 FORMAT('then average value, ENPV is',F20.3
     1/'Std.Dev. ENPV is     ',F26.3)
  853 FORMAT('If set Relative Trigger price at',F15.4)
C**********************************
 1000 ENDDO
C********************************** 
 1100 ENDDO
C**********************************
            
C*********************************
 2000 CONTINUE
      STOP 
      END 
C**********************************