Program PricingDown
!   This code is based on GoofyWelfare.f90       !
!       with paternalistic welfare computed                 !
    
!   Solving Consumer Bankruptcy with Fresh Start    		!
!       with two types of agents - regular and goofy        !
!   Goofy individuals have "wrong" beliefs and don't update !
!       Goffy people believe they have the same shocks      !
!           as regular people, but they don't               !
!   Lenders are not goofy - they do update                  !
!       and keep track of the probability of an individual  !
!       being goofy - that's the type score in this model   !
!   Type score is discetized to keep things tractable       !
!   In bankruptcy, the hospital is paid ahead of the bank   !
!	Note: No shocks in the last period					!
!	No expense shock in the first period				!
!		To change that:									!
!			Fix calculation of EV						!
!			and 1st period in simulations				!
!	Bankruptcy rulled out in the last period			!
!		- computational convenience						!
!	Multiplicative transitory income shocks				!
!		on top of Markovian income process				!
!	Expense shocks proportional to						!
!		Markov component of income process				!
!	Bankruptcy is rulled out immediately following		!
!		a previous bankruptcies.						!
!		Default option is introduced for that period	!


!	Note: Defaulters (after bankruptcy) are not counted as defaulters.	!
!			Their debts do not appear in any of the aggregates.			!

!	Note: Filers of age (lifespan-2) who Default at age (lifespan-1),	!
!			are the only ones allowed to file in last peroid.			!

!	Note: Autonomous consumption is included in measured consumption.	!

!	Warning: Numerical precision during the aggregation may be problematic.		!

Implicit none

Interface
	Real(8) Function Utility(c,cig)
	real(8) c, cig
	end function
end interface


Integer age, avtrans, bin, bmax, count, gridmax, groups, i, id, ii, iii, increason, iscore, isplus, isprime, j
Integer k, kk, kprime, kprimeedu, ks, lifespan, limit                                                    !, limitn, loss
Integer Monte, Monte1, preshock, Pop, PopG, PopR, s, scorebins, scorepoints, sh, shocks, smax, ss, st, sstop, states, trans, transreason
Integer bins(16,5,3), BR(16,5,3,21,6018), BRafter(16,5,3,21,5), BRafteredu(16,5,3,21,5), BRedu(16,5,3,21,6018), DebtRolled(16,5,3,5)
Integer expenseloss(5,5), history(16), kpos(16), PopB(20), scorehis(16), sdummy(1), shockhis(16), stigma(16), transhis(16)

Real(8) asum, asum2, asum3, asum4, aut, autlast, avborint, avborint1, avborintG, avborintR, avdoib, avdoid, avdoibG, avdoibR, avdoidG, avdoidR, avlncon, avlninc
Real(8) baddebt, balance, bankruptcy, bankruptcyG, bankruptcyR, bestluck, betanorm, betalast, borrowers, borrowante, borrowersG, borrowersR, borrowanteG, borrowanteR
Real(8) bprime, burn, burn2, cbar, cigma, costfile, coststig1, coststig2, conalt, countG, countover, countR, countunder
Real(8) debt, debtG, debtR, debtadd, defaults, di, dimaxHS, dimaxLS, disposable, ds, dsmaxHS, dsmaxLS, EV, EVedu, EVpat, gain, garn, gridfine, intinc, lambda, maxborint, maxintgen1
Real(8) medianinc, net, output, outdef, outputG, outputR, outdefG, outdefR, overborrowing, overborrq, overhang, pikplus, pikprime, posterior, prior, profitsG, profitsR, qualscore, qq
Real(8) rb, realizedEV, realizedG, realizedR, recidivist, richfiler, rmaxHS, rmaxLS, rroll, rs, saving, scorebinsize
Real(8) take, tooearly, toolate, varint, varlncon, varlninc, Vcont, Vcontedu, Vcontpat, wage, wannasave, weight1, welfarecum, worstluck, wss				!, Utility
Real(8) A(3), avdebincB(20), avgscoreG(16), avgscoreR(16), avinttotal(100), avintG(100), avintR(100), avlnconage(16), avglnconG(16), avglnconR(16), avlnincage(16)
Real(8) bankrtotal(100), bankrG(100), bankrR(100), beta(16), betacum(16), borpostB(20), borrow(16,5,3,5), borrowG(16,5,3,5), borrowR(16,5,3,5), borrowB(20)
Real(8) borrowbyage(16), borrowbyageG(16), borrowbyageR(16)
Real(8) borrowerposttotal(100), borrowerantetotal(100), borrowerpostG(100), borrowerpostR(100), borroweranteG(100), borroweranteR(100), BRbyreason(5,2,2), BRbyreasonG(5,2,2), BRbyreasonR(5,2,2)
Real(8) cand(6018), candedu(6018), candpat(6018), ccand(6018), Con(16,5,3,21,6018), Con1(16,5,3,21,5), ConB(16,5,3), ConD(16,5,3,5), cons(16), consG(16), consR(16), cutoff(20)
Real(8) debincB(20), debtB(20), debtbyage(16), debtbyageG(16), debtbyageR(16), debtincbortotal(100), debtincdeftotal(100), debtincborG(100), debtincborR(100)       !, debtdist(20,5,5,5), debtincdef(17,5,5,5)
Real(8) debtincdefG(100), debtincdefR(100), debtpostB(20), debttotal(100), debttotalG(100), debttotalR(100), defaulttotal(100)			!, income(3000000), earnings(3000000), wealth(3000000)
Real(8) ES(16), expenconst(5), expenprop(5), expense(5,5), filingsB(20), fracborrow(16,5,3,5), gaintotal(100), hbar(16), incoG(16), incoR(16), interestageG(16), interestageR(16)
Real(8) netprofitG(100), netprofitR(100), newbankruptcy(16), newbankruptcyG(16), newbankruptcyR(16), outputB(20), overborrowtotal(100), overborrqtotal(100), overhangtotal(100)
Real(8) phi(5), probGoofy(5,5), probTrue(5,5), prob0(5), probtGoofy(3), probtTrue(3), probsGoofy(5), probsTrue(5), profitageG(16), profitageR(16)
Real(8) q(16,5,21,6018), realizedwelfare(100), repeatfrac(100), richfrac(100), savebyage(16), savers(100), savetotal(100), totalearly(100), totallate(100), totalover(100), totalunder(100), toolatebyage(16)          !, qtemp(5)
Real(8) V(16,5,3,21,6018), V0(16,5,21), V0alt(16,5,21), V0edu(16,5,21), V0pat(16,5,21), Vaftalt(16,5,3,21,5), Vafter(16,5,3,21,5), Vafteredu(16,5,3,21,5), Vafterpat(16,5,3,21,5)
Real(8) varlnconage(16), varlnconG(16), varlnconR(16), varlnincage(16), Vbaralt(16,5,3,21), Vbaredu(16,5,3,21), VbarFS(16,5,3,21), Vbarpat(16,5,3,21)
Real(8) Vdef(16,5,3,21,5), Vdefalt(16,5,3,21,5), Vdefedu(16,5,3,21,5), Vdefpat(16,5,3,21,5), Vedu(16,5,3,21,6018), Vpat(16,5,3,21,6018)
Real(8) Vprime(16,5,21,6018), Vprimeedu(16,5,21,6018), Vprimepat(16,5,21,6018), welfareG(100), welfareR(100), writeoff(100)   !, welfareexpost(100001)

    
Logical defafter, defedu, flag, goof, run0, runedu, ystig1, ystig2, ystig3
limit=6018
Pop = 100000        !   Population within a single simulation


do i=1,limit
	cand(i)=-1.d96
    candedu(i)=-1.d96
end do
!debtdist=0.0
debtbyage=0.0
debtbyageG=0.0
debtbyageR=0.0
borrow=0.0
borrowG=0.0
borrowR=0.0
fracborrow=0.0
borrowbyage=0.0
borrowbyageG=0.0
borrowbyageR=0.0
savebyage=0.0
toolatebyage=0.0
avgscoreG=0.0
avgscoreR=0.0
PopG=0
PopR=0

!   Parameters  !
open (8,file='inputRevised.txt')
read (8,*) id
read (8,*) lifespan, states, trans, shocks, groups, gridfine, scorebins
read (8,*) cigma, aut, autlast, cbar, rs, rb, rroll, garn, burn, burn2, costfile, coststig1, coststig2, lambda
	!		Elasticity, Autonomous consumption, Autonomous in the last period, Exemption,
	!...	Risk-free savings rate, Risk-free borrowing rate, Roll-over rate
	!...	Fraction of labour income garnished, Fraction of consumption burnt for 2 periods
	!...	Cost of filing, Additive utility stigmas for 2 periods, Share of goofy	    !
read (8,*) qualscore, rmaxLS, rmaxHS, dimaxLS, dimaxHS, dsmaxLS, dsmaxHS
    !       Minimum type score to qualify for loose restrictions, Max borrowing interest rate for low scores,
    !...    Max borrowing interest rate for high scores, Max debt/income for low scores, 
    !...    Max debt/income for high scores, Max debt service for low scores, Max debt service for high scores      !
read (8,*) (hbar(i),i=1,lifespan)				!	Labor endowments			    	!
read (8,*) betanorm, betalast					!	Beta in normal periods, beta for last period	!
read (8,*) (ES(i),i=1,lifespan)					!	Equivalence scales			    	!
read (8,*) (cutoff(i),i=1,(groups-1))			!	Income cutoffs for reporting    	!
read (8,*) (phi(i),i=1,states)					!	Persistent income draws	    		!
read (8,*) (prob0(i),i=1,states)				!	Initial probabilities	    		!
read (8,*) (A(i),i=1,trans)						!	Transitory income draws	    		!
read (8,*) (probtTrue(i),i=1,trans)				!   Realists' Transitory probabilities	!
read (8,*) (probtGoofy(i),i=1,trans)			!	Goofy Transitory probabilities      !
read (8,*) (expenconst(i),i=1,shocks)			!	Expense Shocks - Intercepts     	!
read (8,*) (expenprop(i),i=1,shocks)			!	Expense Shocks - Slopes	    		!
read (8,*) (probsTrue(i),i=1,shocks)            !	Realists' Probabilities of shocks   !
read (8,*) (probsGoofy(i),i=1,shocks)           !	Goofy Probabilities of shocks       !
read (8,*) ((probTrue(i,j),j=1,states),i=1,states)	!	Realists' Transitional probabilities    !
read (8,*) ((probGoofy(i,j),j=1,states),i=1,states)	!	Goofy Transitional probabilities        !
read (8,*) medianinc
    !       Median income                                                               !
close (8)

scorebinsize = 1.0/scorebins
scorepoints=scorebins+1
print *,'The size of a credit score bin is',scorebinsize
! gridfine=500.0

cutoff(groups)=maxval(hbar)*maxval(phi)*maxval(A)+1.0
do age=1,lifespan-1
	do s=1,states
		do st=1,trans
			wage=hbar(age)*phi(s)*A(st)
			i=1
			do while (wage>cutoff(i))
				i=i+1
			end do
			bins(age,s,st)=i
		end do
	end do
end do
do age=1,lifespan-2
	beta(age) = betanorm*(ES(age)/ES(age+1))**(1-cigma)
    if (age>1) then
        betacum(age) = betacum(age-1)*beta(age)
    else
        betacum(age) = beta(age)
    end if
end do
age=lifespan-1
beta(age) = betalast*(ES(age)/ES(age+1))**(1-cigma)
betacum(age) = betacum(age-1)*beta(age)
do i=1,states
	do ii=1,shocks
		expense(i,ii)=expenconst(ii)+expenprop(ii)*phi(i)
	end do
	expense(i,shocks+1)=0.0
end do
probsTrue(shocks+1)=1.0
do sh=1,shocks
	probsTrue(shocks+1)=probsTrue(shocks+1)-probsTrue(sh)
end do
probsGoofy(shocks+1)=1.0
do sh=1,shocks
	probsGoofy(shocks+1)=probsGoofy(shocks+1)-probsGoofy(sh)
end do
avtrans=nint(trans/2.0)

!	Figuring out a grid that fits	!
sstop=0
do while (sstop.eq.0)
	smax=nint(((maxval(phi)-1.0)*(maxval(A)-1.0)*(maxval(hbar)-1.0)*lifespan/2.0)*gridfine)+1
	bmax=nint(smax/2.0)+nint(maxval(expense)*lifespan*gridfine/4.0)+1
!	smax=nint((maxval(phi)*maxval(A)*maxval(hbar)*lifespan/2)*gridfine)+1
!	bmax=smax+nint(maxval(expense)*lifespan*gridfine/2)+1
	gridmax=smax+bmax+1

	print *,'bmax=',bmax
	print *,'gridmax=',gridmax
	if (gridmax>limit) then
		print *,'Grid is too large for arrays to handle. Adjusting grid.'
!		pause
		gridfine=nint(gridfine/2.0)
	else
		sstop=1
	end if
end do
!   Fitting expense shocks to the grid  !
do i=1,states
	do ii=1,shocks
        expenseloss(i,ii)=nint(expense(i,ii)*gridfine)
		expense(i,ii)=expenseloss(i,ii)/gridfine
    end do
    expenseloss(i,shocks+1)=0
end do


print *,'Regular Pr(5,4)=',probTrue(5,4)
print *,'Goofy Pr(5,4)=',probGoofy(5,4)



!	V(age,s,A,score,i) = max [Utility + beta * E[max(V(age+1,s',A',score',i'),Vbar(age+1,s',score',A'))]]	= Value of Repaying
!	Vbar(age,s,A,score) = Utility(burn,garn) + beta * V0(age,s,score)       								= Value of Bankruptcy
!	V0(age,s,score) = E[max(Vafter(age+1,s',A',score',shock'),Vdef(age+1,s',A',score',shock'))]				= Continuation Value during Bankruptcy
!	V0alt(age,s,score) = E[max(Vaftalt(age+1,s',A',score',shock'),Vdefalt(age+1,s',A',score',shock'))]		= Continuation Value during Bankruptcy while saving min possible
!	Vafter(age,s,A,score,shock) = max[Utility(burn2) + beta * Vprime(age,s,score,i')]       				= Value of Repaying after Bankruptcy
!	Vdef(age,s,A,score,shock) = Utility(burn,garn) + beta * Vprime(age,s,score,i'(shock))	        		= Value of Default after Bankruptcy
!	Vprime(age,s,score,i) = E[max(V(age+1,s',A',score',i),Vbar(age+1,s',A',score'))]                        = Continuation Value next period

!	Vedu(age,s,A,score,i) = max [Utility + beta * E[max(Vedu(age+1,s',A',score',i'),Vbaredu(age+1,s',score',A'))]]	= Value of Repaying for self-aware goofy (who will make aware choices from now on)
!	Vbaredu(age,s,A,score) = Utility(burn,garn) + beta * V0edu(age,s,score)       			    					= Value of Bankruptcy for self-aware goofy
!	V0edu(age,s,score) = E[max(Vafteredu(age+1,s',A',score',shock'),Vdefedu(age+1,s',A',score',shock'))]			= Continuation Value during Bankruptcy for self-aware goofy
!	Vafteredu(age,s,A,score,shock) = max[Utility(burn2) + beta * Vprimeedu(age,s,score,i')]         				= Value of Repaying after Bankruptcy for self-aware goofy
!	Vdefedu(age,s,A,score,shock) = Utility(burn,garn) + beta * Vprimeedu(age,s,score,i'(shock))	               		= Value of Default after Bankruptcy for self-aware goofy
!	Vprimeedu(age,s,score,i) = E[max(Vedu(age+1,s',A',score',i),Vbaredu(age+1,s',A',score'))]                       = Continuation Value next period for self-aware goofy

!   All of the Vpat value functions simply keep track of the paternalistic value of making goofy decisions throughout

!   Last Generation !
age=lifespan
do s=1,states
    wage=phi(s)*hbar(age)
	do i=1,gridmax
		net=(i-1-bmax)/gridfine
        if ((wage+net)<0.0) then
	        V(age,s,1:trans,1:scorepoints,i)=-1.d96
	        Vedu(age,s,1:trans,1:scorepoints,i)=-1.d96
	        Vpat(age,s,1:trans,1:scorepoints,i)=-1.d96
		else
			Con(age,s,1:trans,1:scorepoints,i)=wage+net
			V(age,s,1:trans,1:scorepoints,i)=Utility(autlast+wage+net,cigma)
			Vedu(age,s,1:trans,1:scorepoints,i)=Utility(autlast+wage+net,cigma)
			Vpat(age,s,1:trans,1:scorepoints,i)=Utility(autlast+wage+net,cigma)
	    end if
	end do
	ConB(age,s,1:trans)=(wage-max((wage-cbar)*garn,0.0)-costfile)*(1.0-burn)
	if (ConB(age,s,1)<0.0) then
		print *,'Filing costs too high for old guys. State = ',s
!		pause
!		goto 777
	end if
	Con1(age,s,1:trans,1:scorepoints,1:5)=wage*(1.0-burn2)
	V0(age-1,s,1:scorepoints)=Utility(autlast+wage*(1.0-burn2),cigma)-coststig2
	V0edu(age-1,s,1:scorepoints)=Utility(autlast+wage*(1.0-burn2),cigma)-coststig2
	V0pat(age-1,s,1:scorepoints)=Utility(autlast+wage*(1.0-burn2),cigma)-coststig2
    V0alt(age-1,s,1:scorepoints)=Utility(autlast+(wage+1.0/gridfine)*(1.0-burn2),cigma)-coststig2
	!	Note:	There is no possibility of default here:					!
	!			If they went bankrupt last period, their balance is 0		!
	!				and they do not get expense shocks.						!
	print *,V(age,s,1:trans,nint(scorebins/2.0),bmax+1)
	print *,Vedu(age,s,1:trans,nint(scorebins/2.0),bmax+1)
end do
!pause

!   Other Generations   !
do age=(lifespan-1),1,-1
    do s=1,states

		!   Pricing the debt (using score-weighted p) and Calculating continuation values  !
		do k=1,gridmax
            balance=(k-1-bmax)/gridfine
            do iscore=1,scorepoints
                prior=(iscore-1.0)*scorebinsize                
    			if (age<lifespan-1) then
	    			Vcont=0.0
	    			Vcontedu=0.0
	    			Vcontpat=0.0
		    		qq=0.0
                    asum=0.0
                    asum2=0.0
                    asum3=0.0
			    	do ss=1,states
				    	do st=1,trans
    						wss=phi(ss)*A(st)*hbar(age+1)
	    					do sh=1,shocks+1
!		    					loss=nint(expense(ss,sh)*gridfine)
			    				ks=max(k-expenseloss(ss,sh),1)
                                posterior = (prior*probTrue(s,ss)*probtTrue(st)*ProbsTrue(sh))/(prior*probTrue(s,ss)*probtTrue(st)*ProbsTrue(sh) + (1-prior)*probGoofy(s,ss)*probtGoofy(st)*ProbsGoofy(sh))
                                isprime = int(posterior/scorebinsize)+1
                                isplus = min(isprime+1,scorepoints)
                                weight1 = (isprime*scorebinsize - posterior)/scorebinsize
                                !       Since posterior is between score points, the score is randomized appropriately:
                                !           It is assigned ISprime with probability weight1, and ISplus with remaining probability
				    			if (V(age+1,ss,st,isprime,ks)<VbarFS(age+1,ss,st,isprime)) then
					    			Vcont=Vcont+weight1*probTrue(s,ss)*probtTrue(st)*probsTrue(sh)*VbarFS(age+1,ss,st,isprime)
					    			Vcontpat=Vcontpat+weight1*probGoofy(s,ss)*probtGoofy(st)*probsGoofy(sh)*Vbarpat(age+1,ss,st,isprime)
					    			asum=asum+weight1*probTrue(s,ss)*probtTrue(st)*probsTrue(sh)
						    		if (balance<0.0) then
							    		qq=qq+max((wss-expense(ss,sh)-cbar)*garn,0.0)*weight1*(prior*probTrue(s,ss)*probtTrue(st)*ProbsTrue(sh) + (1-prior)*probGoofy(s,ss)*probtGoofy(st)*ProbsGoofy(sh))/(-balance)
							    		asum2=asum2+weight1*(prior*probTrue(s,ss)*probtTrue(st)*ProbsTrue(sh) + (1-prior)*probGoofy(s,ss)*probtGoofy(st)*ProbsGoofy(sh))
								    else
									    qq=qq+weight1*(prior*probTrue(s,ss)*probtTrue(st)*ProbsTrue(sh) + (1-prior)*probGoofy(s,ss)*probtGoofy(st)*ProbsGoofy(sh))
							    		asum2=asum2+weight1*(prior*probTrue(s,ss)*probtTrue(st)*ProbsTrue(sh) + (1-prior)*probGoofy(s,ss)*probtGoofy(st)*ProbsGoofy(sh))
    								end if
	    						else
		    						Vcont=Vcont+weight1*probTrue(s,ss)*probtTrue(st)*probsTrue(sh)*V(age+1,ss,st,isprime,ks)
		    						Vcontpat=Vcontpat+weight1*probGoofy(s,ss)*probtGoofy(st)*probsGoofy(sh)*Vpat(age+1,ss,st,isprime,ks)
					    			asum=asum+weight1*probTrue(s,ss)*probtTrue(st)*probsTrue(sh)
			    					qq=qq+weight1*(prior*probTrue(s,ss)*probtTrue(st)*ProbsTrue(sh) + (1-prior)*probGoofy(s,ss)*probtGoofy(st)*ProbsGoofy(sh))
                                    asum2=asum2+weight1*(prior*probTrue(s,ss)*probtTrue(st)*ProbsTrue(sh) + (1-prior)*probGoofy(s,ss)*probtGoofy(st)*ProbsGoofy(sh))
				    			end if
				    			if (Vedu(age+1,ss,st,isprime,ks)<Vbaredu(age+1,ss,st,isprime)) then
					    			Vcontedu=Vcontedu+weight1*probGoofy(s,ss)*probtGoofy(st)*probsGoofy(sh)*Vbaredu(age+1,ss,st,isprime)
	    						else
		    						Vcontedu=Vcontedu+weight1*probGoofy(s,ss)*probtGoofy(st)*probsGoofy(sh)*Vedu(age+1,ss,st,isprime,ks)
				    			end if
				    			if (V(age+1,ss,st,isplus,ks)<VbarFS(age+1,ss,st,isplus)) then
					    			Vcont=Vcont+(1-weight1)*probTrue(s,ss)*probtTrue(st)*probsTrue(sh)*VbarFS(age+1,ss,st,isplus)
					    			Vcontpat=Vcontpat+(1-weight1)*probGoofy(s,ss)*probtGoofy(st)*probsGoofy(sh)*Vbarpat(age+1,ss,st,isplus)
					    			asum=asum+(1-weight1)*probTrue(s,ss)*probtTrue(st)*probsTrue(sh)
						    		if (balance<0.0) then
							    		qq=qq+max((wss-expense(ss,sh)-cbar)*garn,0.0)*(1-weight1)*(prior*probTrue(s,ss)*probtTrue(st)*ProbsTrue(sh) + (1-prior)*probGoofy(s,ss)*probtGoofy(st)*ProbsGoofy(sh))/(-balance)
									    asum2=asum2+(1-weight1)*(prior*probTrue(s,ss)*probtTrue(st)*ProbsTrue(sh) + (1-prior)*probGoofy(s,ss)*probtGoofy(st)*ProbsGoofy(sh))
								    else
									    qq=qq+(1-weight1)*(prior*probTrue(s,ss)*probtTrue(st)*ProbsTrue(sh) + (1-prior)*probGoofy(s,ss)*probtGoofy(st)*ProbsGoofy(sh))
									    asum2=asum2+(1-weight1)*(prior*probTrue(s,ss)*probtTrue(st)*ProbsTrue(sh) + (1-prior)*probGoofy(s,ss)*probtGoofy(st)*ProbsGoofy(sh))
    								end if
	    						else
		    						Vcont=Vcont+(1-weight1)*probTrue(s,ss)*probtTrue(st)*probsTrue(sh)*V(age+1,ss,st,isplus,ks)
		    						Vcontpat=Vcontpat+(1-weight1)*probGoofy(s,ss)*probtGoofy(st)*probsGoofy(sh)*Vpat(age+1,ss,st,isplus,ks)
					    			asum=asum+(1-weight1)*probTrue(s,ss)*probtTrue(st)*probsTrue(sh)
			    					qq=qq+(1-weight1)*(prior*probTrue(s,ss)*probtTrue(st)*ProbsTrue(sh) + (1-prior)*probGoofy(s,ss)*probtGoofy(st)*ProbsGoofy(sh))
                                    asum2=asum2+(1-weight1)*(prior*probTrue(s,ss)*probtTrue(st)*ProbsTrue(sh) + (1-prior)*probGoofy(s,ss)*probtGoofy(st)*ProbsGoofy(sh))
				    			end if
				    			if (Vedu(age+1,ss,st,isplus,ks)<Vbaredu(age+1,ss,st,isplus)) then
					    			Vcontedu=Vcontedu+(1-weight1)*probGoofy(s,ss)*probtGoofy(st)*probsGoofy(sh)*Vbaredu(age+1,ss,st,isplus)
	    						else
		    						Vcontedu=Vcontedu+(1-weight1)*probGoofy(s,ss)*probtGoofy(st)*probsGoofy(sh)*Vedu(age+1,ss,st,isplus,ks)
				    			end if
					    	end do
    		            end do
	    			end do
		    	else
			    	wss=phi(s)*hbar(age+1)
				    if (wss+balance<0.0) then
    					Vcont=-1.d96
    					Vcontpat=-1.d96
                        asum=1.0
    					Vcontedu=-1.d96
	    				qq=0.0
                        asum2=1.0
		    		else
			    		Vcont=V(age+1,s,1,1,k)
			    		Vcontpat=Vpat(age+1,s,1,1,k)
                        asum=1.0
			    		Vcontedu=Vedu(age+1,s,1,1,k)
				    	qq=1.0
                        asum2=1.0
    				end if
	    		end if

    			if (balance<0.0) then
                    di = -balance*qq/(phi(s)*hbar(age)*(1.0+rb)) 
                    ds = -(1-qq/(1.0+rb))*balance/(phi(s)*hbar(age))
                    if (prior<qualscore) then
    	    			if (qq<(1.0+rb)/(1.0+rmaxLS)) then		!	Interest rate ceiling binds	!
	    	    			q(age,s,iscore,k)=0.0
                        elseif (di>dimaxLS) then
                            q(age,s,iscore,k)=0.0
                        elseif (ds>dsmaxLS) then
                            q(age,s,iscore,k)=0.0
		    	    	else
			    	    	q(age,s,iscore,k)=qq/(1.0+rb)
                        end if
                    else
    	    			if (qq<(1.0+rb)/(1.0+rmaxHS)) then		!	Interest rate ceiling binds	!
	    	    			q(age,s,iscore,k)=0.0
                        elseif (di>dimaxHS) then
                            q(age,s,iscore,k)=0.0
                        elseif (ds>dsmaxHS) then
                            q(age,s,iscore,k)=0.0
		    	    	else
			    	    	q(age,s,iscore,k)=qq/(1.0+rb)
                        end if
                    end if
	    		else
		    		q(age,s,iscore,k)=1.0/(1.0+rs)
			    end if
    			Vprime(age,s,iscore,k)=Vcont
    			Vprimeedu(age,s,iscore,k)=Vcontedu
    			Vprimepat(age,s,iscore,k)=Vcontpat
                if (abs(asum-1.0)>0.0000001) then
                    print *,'Probabilities do not add up. Value calculation adds to',asum,'in age',age,'state',s,'type-score',iscore,'future assets',k
                    pause
                end if
                if (abs(asum2-1.0)>0.0000001) then
                    print *,'Probabilities do not add up. Price calculation adds to',asum2,'in age',age,'state',s,'type-score',iscore,'future assets',k
                    pause
                end if
                
            end do
		end do


		print *,'State',s


		do st=1,trans
		    wage=phi(s)*A(st)*hbar(age)

!			print *,'Transitory',st

	        do iscore=1,scorepoints
                do i=1,gridmax
    		        net=(i-1-bmax)/gridfine
    			    do k=1,gridmax
	    			    balance=(k-1-bmax)/gridfine
		    			disposable=net-q(age,s,iscore,k)*balance
	                    if ((wage+disposable)<=0.0) then
		                    cand(k)=-1.d96
			            else
    						ccand(k)=wage+disposable
	    					cand(k)=Utility(aut+wage+disposable,cigma)
	                    end if
		                candedu(k)=cand(k)+beta(age)*Vprimeedu(age,s,iscore,k)
                        candpat(k)=cand(k)+beta(age)*Vprimepat(age,s,iscore,k)
		                cand(k)=cand(k)+beta(age)*Vprime(age,s,iscore,k)
    			    end do

    	            !   Maximize    !
	    	        V(age,s,st,iscore,i)=maxval(cand)
    				sdummy=maxloc(cand)
	    			BR(age,s,st,iscore,i)=sdummy(1)
    				Con(age,s,st,iscore,i)=ccand(sdummy(1))
                    Vpat(age,s,st,iscore,i) = candpat(sdummy(1))

	    	        Vedu(age,s,st,iscore,i)=maxval(candedu)
    				sdummy=maxloc(candedu)
	    			BRedu(age,s,st,iscore,i)=sdummy(1)
    	        end do

    			if (age>1) then
	    			do sh=1,shocks+1
		    			net=-expense(s,sh)

    					!	Value of Repaying after Bankruptcy	!
	    	            do k=1,gridmax
    			            balance=(k-1-bmax)/gridfine
	    			        disposable=net-q(age,s,iscore,k)*balance
		    			    if ((wage+disposable)<=0.0) then
			    			    cand(k)=-1.d96
		                    else
    							ccand(k)=(wage+disposable)*(1-burn2)
	    						cand(k)=Utility(aut+(wage+disposable)*(1.0-burn2),cigma)
		    			    end if
			    			candedu(k)=cand(k)+beta(age)*Vprimeedu(age,s,iscore,k)
			    			candpat(k)=cand(k)+beta(age)*Vprimepat(age,s,iscore,k)
			    			cand(k)=cand(k)+beta(age)*Vprime(age,s,iscore,k)
    					end do

    		            !   Maximize    !
	    		        Vafter(age,s,st,iscore,sh)=maxval(cand)-coststig2
    					sdummy=maxloc(cand)
	    				BRafter(age,s,st,iscore,sh)=sdummy(1)
    					Con1(age,s,st,iscore,sh)=ccand(sdummy(1))
                        Vafterpat(age,s,st,iscore,sh)=candpat(sdummy(1))

	    		        Vafteredu(age,s,st,iscore,sh)=maxval(candedu)-coststig2
    					sdummy=maxloc(candedu)
	    				BRafteredu(age,s,st,iscore,sh)=sdummy(1)


    					!	Value of Repaying after Bankruptcy with Saving 1	!
	    	            do k=1,gridmax
		    	            balance=(k-1-bmax)/gridfine
			    	        disposable=net+1.0/gridfine - q(age,s,iscore,k)*balance
    					    if ((wage+disposable)<=0.0) then
	    					    cand(k)=-1.d96
		                    else
			    				ccand(k)=(wage+disposable)*(1-burn2)
				    			cand(k)=Utility(aut+(wage+disposable)*(1.0-burn2),cigma)
    					    end if
	    					cand(k)=cand(k)+beta(age)*Vprime(age,s,iscore,k)
    					end do
	    	            !   Maximize    !
		    	        Vaftalt(age,s,st,iscore,sh)=maxval(cand)-coststig2


    					!	Value of Default after Bankruptcy	!
	    				take=max((wage-cbar)*garn,0.0)
		    			bprime=min((net+take)*(1+rroll),0.0)
			    		DebtRolled(age,s,st,sh)=bmax+1+nint(bprime*gridfine)
    					ConD(age,s,st,sh)=(wage-take)*(1.0-burn)
    					if (age<lifespan-1) then
	    					Vdef(age,s,st,iscore,sh)=Utility(aut+ConD(age,s,st,sh),cigma)-coststig1+beta(age)*Vprime(age,s,iscore,DebtRolled(age,s,st,sh))
	    					Vdefpat(age,s,st,iscore,sh)=Utility(aut+ConD(age,s,st,sh),cigma)-coststig1+beta(age)*Vprimepat(age,s,iscore,DebtRolled(age,s,st,sh))
	    					Vdefedu(age,s,st,iscore,sh)=Utility(aut+ConD(age,s,st,sh),cigma)-coststig1+beta(age)*Vprimeedu(age,s,iscore,DebtRolled(age,s,st,sh))
    					else
	    					Vdef(age,s,st,iscore,sh)=Utility(aut+ConD(age,s,st,sh),cigma)-coststig1+beta(age)*(Utility(autlast+ConB(age+1,s,1),cigma)-coststig1)
	    					Vdefpat(age,s,st,iscore,sh)=Utility(aut+ConD(age,s,st,sh),cigma)-coststig1+beta(age)*(Utility(autlast+ConB(age+1,s,1),cigma)-coststig1)
	    					Vdefedu(age,s,st,iscore,sh)=Utility(aut+ConD(age,s,st,sh),cigma)-coststig1+beta(age)*(Utility(autlast+ConB(age+1,s,1),cigma)-coststig1)
                        !       We are forcing bankruptcy in retirement period if there was bankruptcy at (J-2) AND expense shock at (J-1) was defaulted on
    					end if

    					!	Value of Default after Bankruptcy with Saving 1	!
	    				take=max((wage-cbar)*garn,0.0)
		    			bprime=min((net+1.0/gridfine+take)*(1+rroll),0.0)
			    		conalt=(wage-take)*(1.0-burn)
				    	if (age<lifespan-1) then
					    	Vdefalt(age,s,st,iscore,sh)=Utility(aut+conalt,cigma)-coststig1+beta(age)*Vprime(age,s,iscore,bmax+1+nint(bprime*gridfine))
    					else
	    					Vdefalt(age,s,st,iscore,sh)=Utility(aut+conalt,cigma)-coststig1+beta(age)*(Utility(autlast+ConB(age+1,s,1),cigma)-coststig1)
    					end if
    			    end do
                end if
            end do            

	        !   Value of Bankruptcy   !
			ConB(age,s,st)=(wage-max((wage-cbar)*garn,0.0)-costfile)*(1.0-burn)
			if (ConB(age,s,st)<0.0) then
				print *,'Filing costs too high. Age, persist, trans = ',age,s,st
				pause
				goto 777
            end if
            do iscore=1,scorepoints
			    VbarFS(age,s,st,iscore)=Utility(aut+ConB(age,s,st),cigma)-coststig1+beta(age)*V0(age,s,iscore)
			    Vbarpat(age,s,st,iscore)=Utility(aut+ConB(age,s,st),cigma)-coststig1+beta(age)*V0pat(age,s,iscore)
			    Vbaredu(age,s,st,iscore)=Utility(aut+ConB(age,s,st),cigma)-coststig1+beta(age)*V0edu(age,s,iscore)
            end do            

	        !   Value of Bankruptcy with Saving 1	!
			conalt=(wage-max((wage-cbar)*garn,0.0)-costfile-1.0/(gridfine*(1.0+rs)))*(1.0-burn)
			if (conalt<0.0) then
				print *,'Filing costs almost too high. Age, persist, trans = ',age,s,st
				pause
				goto 777
			end if
            do iscore=1,scorepoints
    			Vbaralt(age,s,st,iscore)=Utility(aut+conalt,cigma)-coststig1+beta(age)*V0alt(age,s,iscore)
            end do
	    end do
	end do

    if (age>1) then
        do s=1,states
            do iscore=1,scorepoints
                prior=(iscore-1.0)*scorebinsize                
                asum=0.0
	    		asum2=0.0
                asum3=0.0
                asum4=0.0
                do ss=1,states
			    	do st=1,trans
				    	do sh=1,shocks+1
                            posterior = (prior*probTrue(s,ss)*probtTrue(st)*ProbsTrue(sh))/(prior*probTrue(s,ss)*probtTrue(st)*ProbsTrue(sh) + (1-prior)*probGoofy(s,ss)*probtGoofy(st)*ProbsGoofy(sh))
                            isprime = int(posterior/scorebinsize)+1
                            isplus = min(isprime+1,scorepoints)
                            weight1 = (isprime*scorebinsize - posterior)/scorebinsize
                            if (weight1<-0.01) then
                                print *,'Negative weight=',weight1
                                pause
                            end if
                            if (Vdef(age,ss,st,isprime,sh)>Vafter(age,ss,st,isprime,sh)) then
                                pikprime = Vdefpat(age,ss,st,isprime,sh)
                            else
                                pikprime = Vafterpat(age,ss,st,isprime,sh)
                            end if
                            if (Vdef(age,ss,st,isplus,sh)>Vafter(age,ss,st,isplus,sh)) then
                                pikplus = Vdefpat(age,ss,st,isplus,sh)
                            else
                                pikplus = Vafterpat(age,ss,st,isplus,sh)
                            end if                            
asum=asum+weight1*probTrue(s,ss)*probtTrue(st)*probsTrue(sh)*max(Vafter(age,ss,st,isprime,sh),Vdef(age,ss,st,isprime,sh))+(1-weight1)*probTrue(s,ss)*probtTrue(st)*probsTrue(sh)*max(Vafter(age,ss,st,isplus,sh),Vdef(age,ss,st,isplus,sh))
asum2=asum2+weight1*probTrue(s,ss)*probtTrue(st)*probsTrue(sh)*max(Vaftalt(age,ss,st,isprime,sh),Vdefalt(age,ss,st,isprime,sh))+(1-weight1)*probTrue(s,ss)*probtTrue(st)*probsTrue(sh)*max(Vaftalt(age,ss,st,isplus,sh),Vdefalt(age,ss,st,isplus,sh))
asum3=asum3+weight1*probGoofy(s,ss)*probtGoofy(st)*probsGoofy(sh)*max(Vafteredu(age,ss,st,isprime,sh),Vdefedu(age,ss,st,isprime,sh))+(1-weight1)*probGoofy(s,ss)*probtGoofy(st)*probsGoofy(sh)*max(Vafteredu(age,ss,st,isplus,sh),Vdefedu(age,ss,st,isplus,sh))
asum4=asum4+weight1*probGoofy(s,ss)*probtGoofy(st)*probsGoofy(sh)*pikprime+(1-weight1)*probGoofy(s,ss)*probtGoofy(st)*probsGoofy(sh)*pikplus
	    				end do
		            end do
			    end do
                V0(age-1,s,iscore)=asum
                V0alt(age-1,s,iscore)=asum2
                V0edu(age-1,s,iscore)=asum3
                V0pat(age-1,s,iscore)=asum4
            end do
        end do
    end if


	print *,'Age=',age
	do s=1,states
		print *,V(age,s,1:trans,nint(scorebins/2.0),bmax+1)
		print *,Vedu(age,s,1:trans,nint(scorebins/2.0),bmax+1)
	end do
	do s=1,states
		print *,BR(age,s,1:trans,nint(scorebins/2.0),bmax+1)
	end do
    !pause
end do

EV=0.0
EVpat=0.0
EVedu=0.0
do s=1,states
	do st=1,trans
        posterior = ((1-lambda)*probtTrue(st))/((1-lambda)*probtTrue(st) + lambda*probtGoofy(st))
        iscore = min(int(posterior/scorebinsize)+1,scorepoints)
        isplus = min(iscore+1,scorepoints)
        weight1 = (iscore*scorebinsize - posterior)/scorebinsize
	    EV=EV+weight1*prob0(s)*probtTrue(st)*V(1,s,st,iscore,bmax+1)+(1-weight1)*prob0(s)*probtTrue(st)*V(1,s,st,isplus,bmax+1)
	    EVpat=EVpat+weight1*prob0(s)*probtGoofy(st)*Vpat(1,s,st,iscore,bmax+1)+(1-weight1)*prob0(s)*probtGoofy(st)*Vpat(1,s,st,isplus,bmax+1)
	    EVedu=EVedu+weight1*prob0(s)*probtGoofy(st)*Vedu(1,s,st,iscore,bmax+1)+(1-weight1)*prob0(s)*probtGoofy(st)*Vedu(1,s,st,isplus,bmax+1)
	end do
end	do
print *,'Realist EV=',EV
print *,'Paternalistic Goofy EV=',EVpat
print *,'Educated Goofy EV=',EVedu
!       EVedu is not the paternalistic welfare. EVedu computed the value of goofy people making correct AWARE choices

999 format (F8.4,3F8.5)
open (21,file='PricesYoung.txt')
write (21,*) 'Input ID#',id
write (21,*) 'Age=1'
ii=nint(scorebins/2.0)
write (21,*) 'LowScore is', 1.0/scorebins
write (21,*) 'MedianScore is', (ii-1.0)/scorebins
write (21,*) 'HighScore is', (scorebins-1.0)/scorebins
write (21,*) 'Balance ','PricesLowScore ','PricesMedianScore ','PricesHighScore '
do i=1,(bmax+3)
	net=(i-1-bmax)/gridfine
	write (21,999) net,q(1,1,2,i),q(1,1,ii,i),q(1,1,scorebins,i)
end do
close (21)
open (24,file='PricesMiddleAge.txt')
write (24,*) 'Input ID#',id
age=nint(lifespan*0.35)
write (24,*) 'Age=',age
write (24,*) 'LowScore is', 1.0/scorebins
write (24,*) 'MedianScore is', (ii-1.0)/scorebins
write (24,*) 'HighScore is', (scorebins-1.0)/scorebins
write (24,*) 'Balance ','PricesLowScore ','PricesMedianScore ','PricesHighScore '
do i=1,(bmax+3)
	net=(i-1-bmax)/gridfine
	write (24,999) net,q(age,1,2,i),q(age,1,ii,i),q(age,1,scorebins,i)
end do
close (24)
open (25,file='PricesOld.txt')
write (25,*) 'Input ID#',id
age=nint(lifespan*0.85)
write (25,*) 'Age=',age
write (25,*) 'LowScore is', 1.0/scorebins
write (25,*) 'MedianScore is', (ii-1.0)/scorebins
write (25,*) 'HighScore is', (scorebins-1.0)/scorebins
write (25,*) 'Balance ','PricesLowScore ','PricesMedianScore ','PricesHighScore '
do i=1,(bmax+3)
	net=(i-1-bmax)/gridfine
	write (25,999) net,q(age,1,2,i),q(age,1,ii,i),q(age,1,scorebins,i)
end do
close (25)
!
!998 format (F8.4,15F8.5)
!open (22,file='qFSage.txt')
!write (22,*) 'Input ID#',id
!write (22,*) 'EV=',EV
!write (22,*) 'Prices_for_individuals_with_persistent_state_3_by_age'
!write (22,*) 'Balance ','Prices_by_age'
!do i=1,(bmax+3)
!	net=(i-1-bmax)/gridfine
!	write (22,998) net,q(1:(lifespan-1),3,i)
!end do
!close (22)



!   			Simulations								    !
!   	Pop = Population per generation in a simulation 	!
!   	Monte = Number of simulations					    !

Monte1 = 15
Monte = 30
iii = 8657658		!   	Initial seed for RNG    	!

!998 format (2I3,3F8.5)
!open (26,file='DebtDumpGoofy.txt')
!write (26,*) 'Input ID#',id
!write (26,*) 'Age ','ScoreBin ','AmountBorrowed ','Price ','FullInfoPrice'
!open (27,file='DebtDumpRealist.txt')
!write (27,*) 'Input ID#',id
!write (27,*) 'Age ','ScoreBin ','AmountBorrowed ','Price ','FullInfoPrice'

!count=0
do age=1,lifespan
	newbankruptcy(age)=0.0
	newbankruptcyG(age)=0.0
	newbankruptcyR(age)=0.0
	cons(age)=0.0
	avlnconage(age)=0.0
	avglnconG(age)=0.0    
	avglnconR(age)=0.0    
	avlnincage(age)=0.0
end do
BRbyreason(1:5,1:2,1:2)=0.0
BRbyreasonG(1:5,1:2,1:2)=0.0
BRbyreasonR(1:5,1:2,1:2)=0.0
outdef=0.0
output=0.0
outdefG=0.0
outputG=0.0
outdefR=0.0
outputR=0.0
incoG=0.0
incoR=0.0
interestageG=0.0
interestageR=0.0
profitageG=0.0
profitageR=0.0

PopB=0
outputB=0.0
debtB=0.0
borrowB=0.0
borpostB=0.0
debtpostB=0.0
avdebincB=0.0
filingsB=0.0

bestluck = -1.d96
worstluck = 1.d96

do j=1,Monte1
	debt=0.0
	debtG=0.0
	debtR=0.0
	saving=0.0
	avborint=0.0
	avborintG=0.0
	avborintR=0.0
	borrowante=0.0
	borrowanteG=0.0
	borrowanteR=0.0
	bankruptcy=0.0
    recidivist=0.0
    richfiler=0.0
	bankruptcyG=0.0
	bankruptcyR=0.0
	wannasave=0.0
	defaults=0.0
	borrowers=0.0
	borrowersG=0.0
	borrowersR=0.0
	avdoib=0.0
	avdoibG=0.0
	avdoibR=0.0
	avdoid=0.0
	avdoidG=0.0
	avdoidR=0.0
	baddebt=0.0
    realizedEV=0.0
    realizedG=0.0
    realizedR=0.0
    overborrowing=0.0
    overborrq=0.0
    overhang=0.0
    tooearly=0.0
    toolate=0.0
    gain=0.0
    profitsG=0.0
    profitsR=0.0
    countG = 0.0
    countR = 0.0
    countover = 0.0
    countunder = 0.0
    print *,'Monte=',j


	do i=1,Pop
        welfarecum=0.0
        flag=.FALSE.
        age=1
		wss=ran(iii)
		!	Assigning type	!
        goof=(wss<lambda)
        if (goof) then
            PopG=PopG+1
            countG=countG+1.0
        else
            PopR=PopR+1
            countR=countR+1.0
        endif     

        wss=ran(iii)
		!	Assigning Markov productivity	!
		sstop=0
		asum=0.0
		s=1
		do while (sstop.eq.0)
			if (wss < asum + prob0(s)) then
				history(1)=s
				sstop=1
				if (s>states) then
					print *,'Run #',j,' Markov =',s,' Random =',wss
					history(1)=states
				end if
			else
				asum = asum + prob0(s)
				s=s+1
			end if
		end do

		wss=ran(iii)
		!	Assigning transitory productivity	!
		sstop=0
		asum=0.0
		st=1
        if (goof) then
    		do while (sstop.eq.0)
	    		if (wss < asum + probtGoofy(st)) then
    				transhis(1)=st
	    			sstop=1
		    		if (st>trans) then
			    		print *,'Run #',j,' Transitory =',st,' Random =',wss
				    	transhis(1)=trans
    				end if
	    		else
		    		asum = asum + probtGoofy(st)
    				st=st+1
	    		end if
            end do
        else   
		    do while (sstop.eq.0)
			    if (wss < asum + probtTrue(st)) then
				    transhis(1)=st
    				sstop=1
	    			if (st>trans) then
		    			print *,'Run #',j,' Transitory =',st,' Random =',wss
			    		transhis(1)=trans
				    end if
    			else
	    			asum = asum + probtTrue(st)
		    		st=st+1
    			end if
            end do
        end if

        posterior = ((1-lambda)*probtTrue(st))/((1-lambda)*probtTrue(st) + lambda*probtGoofy(st))
        isprime = min(int(posterior/scorebinsize)+1,scorepoints)
        weight1 = (isprime*scorebinsize - posterior)/scorebinsize
        wss=ran(iii)
        if (wss<weight1) then
            scorehis(1) = isprime
        else
            scorehis(1) = min(isprime+1,scorepoints)
        end if
        
        if (goof) then
            avgscoreG(1)=avgscoreG(1)+scorehis(1)-1
        else
            avgscoreR(1)=avgscoreR(1)+scorehis(1)-1
        end if

		shockhis(1)=shocks+1
		kpos(1)=bmax+1        

		wage=phi(history(1))*A(transhis(1))*hbar(1)
		bin=bins(1,history(1),transhis(1))
		PopB(bin)=PopB(bin)+1
		outputB(bin)=outputB(bin)+wage

		if (wage>0.0) then
			output=output+wage
            if (goof) then
                outputG=outputG+wage
                incoG(1)=incoG(1)+wage
            else
                outputR=outputR+wage
                incoR(1)=incoR(1)+wage
            endif
			avlnincage(1)=avlnincage(1)+log(wage)
		else
			avlnincage(1)=avlnincage(1)-1.0d80
		end if

        stigma(1)=0
        kprime=BR(1,history(1),transhis(1),scorehis(1),kpos(1))
        cons(1)=cons(1)+Con(1,history(1),transhis(1),scorehis(1),kpos(1))+aut
        avlnconage(1)=avlnconage(1)+log(Con(1,history(1),transhis(1),scorehis(1),kpos(1))+aut)
        if (goof) then
            avglnconG(1)=avglnconG(1)+log(Con(1,history(1),transhis(1),scorehis(1),kpos(1))+aut)
        else
            avglnconR(1)=avglnconR(1)+log(Con(1,history(1),transhis(1),scorehis(1),kpos(1))+aut)
        end if
        welfarecum=welfarecum + Utility(aut+Con(1,history(1),transhis(1),scorehis(1),kpos(1)),cigma)
        
        if (goof) then
            kprimeedu=BRedu(1,history(1),transhis(1),scorehis(1),kpos(1))
            gain = gain + Vedu(1,history(1),transhis(1),scorehis(1),kpos(1)) - (V(1,history(1),transhis(1),scorehis(1),kpos(1)) - beta(1)*Vprime(1,history(1),scorehis(1),kprime) + beta(1)*Vprimeedu(1,history(1),scorehis(1),kprime))
            if (kprime<min(kprimeedu,bmax+1)) then
                overborrowing = overborrowing + (min(kprimeedu,bmax+1)-kprime)/gridfine
                overborrq = overborrq + q(1,history(1),scorehis(1),kprime)*(min(kprimeedu,bmax+1)-kprime)/gridfine
                countover=countover+1.0
            elseif (kprimeedu<min(kprime,bmax+1)) then
                countunder=countunder+1.0
            end if
        end if

		if (kprime<bmax+1) then
			debtadd=-(q(1,history(1),scorehis(1),kprime)*(kprime-1-bmax)/gridfine)
			if (debtadd<0) then
				print *,'Debt<0'
				print *,'q=',q(1,history(1),scorehis(1),kprime)
				!pause
			else
				debt=debt+debtadd
				debtB(bin)=debtB(bin)+debtadd
				borrowB(bin)=borrowB(bin)+1.0
				debtbyage(1)=debtbyage(1)+debtadd
				borrowante=borrowante+1.0
				borrow(1,history(1),transhis(1),shockhis(1))=borrow(1,history(1),transhis(1),shockhis(1)) + 1.0
				avborint=avborint-(1-q(1,history(1),scorehis(1),kprime))*(kprime-bmax-1)/gridfine
                if (goof) then
    				debtG=debtG+debtadd
                    debtbyageG(1)=debtbyageG(1)+debtadd
    				borrowanteG=borrowanteG+1.0
    				borrowG(1,history(1),transhis(1),shockhis(1))=borrowG(1,history(1),transhis(1),shockhis(1)) + 1.0
    				avborintG=avborintG-(1-q(1,history(1),scorehis(1),kprime))*(kprime-bmax-1)/gridfine
    				interestageG(1)=interestageG(1)-(1-q(1,history(1),scorehis(1),kprime))*(kprime-bmax-1)/gridfine
                    profitsG=profitsG-debtadd
                    profitageG(1)=profitageG(1)-debtadd
!                    write (26,998) age,scorehis(1),debtadd,q(1,history(1),scorehis(1),kprime),q(1,history(1),1,kprime)
                else
    				debtR=debtR+debtadd
                    debtbyageR(1)=debtbyageR(1)+debtadd
    				borrowanteR=borrowanteR+1.0
    				borrowR(1,history(1),transhis(1),shockhis(1))=borrowR(1,history(1),transhis(1),shockhis(1)) + 1.0
    				avborintR=avborintR-(1-q(1,history(1),scorehis(1),kprime))*(kprime-bmax-1)/gridfine
    				interestageR(1)=interestageR(1)-(1-q(1,history(1),scorehis(1),kprime))*(kprime-bmax-1)/gridfine
                    profitsR=profitsR-debtadd
                    profitageR(1)=profitageR(1)-debtadd
!                    write (27,998) age,scorehis(1),debtadd,q(1,history(1),scorehis(1),kprime),q(1,history(1),scorepoints,kprime)
                end if
			end if
		else
			debtadd=((kprime-1-bmax)/gridfine)/(1+rs)
			saving=saving+debtadd
			savebyage(1)=savebyage(1)+debtadd
		end if


		do age=2,lifespan-1
            if (goof) then
    			wss=ran(iii)
	    		!	Assigning Markov productivity	!
		    	sstop=0
			    asum=0.0
    			s=1
	    		do while (sstop.eq.0)
		    		if (wss < asum + probGoofy(history(age-1),s)) then
			    		history(age)=s
				    	sstop=1
					    if (s>states) then
    						print *,'Run #',j,' Person #',i,' Age =',age,' Markov =',s,' Random =',wss
	    					history(age)=states
		    			end if
    				else
	    				asum = asum + probGoofy(history(age-1),s)
		    			s=s+1
			    	end if
    			end do

    			wss=ran(iii)
	    		!	Assigning transitory productivity	!
		    	sstop=0
			    asum=0.0
    			st=1
	    		do while (sstop.eq.0)
		    		if (wss < asum + probtGoofy(st)) then
			    		transhis(age)=st
				    	sstop=1
    					if (st>trans) then
	    					print *,'Run #',j,' Person #',i,' Age =',age,' Transitory =',st,' Random =',wss
		    				transhis(age)=trans
			    		end if
				    else
					    asum = asum + probtGoofy(st)
    					st=st+1
	    			end if
		    	end do

    			wss=ran(iii)
	    		!	Assigning expense shock	!
		    	sstop=0
			    asum=0.0
    			sh=1
	    		do while (sstop.eq.0)
		    		if (wss < asum + probsGoofy(sh)) then
			    		shockhis(age)=sh
				    	sstop=1
					    if (sh>shocks+1) then
    						print *,'Run #',j,' Person #',i,' Age =',age,' Expense =',sh,' Random =',wss
	    					shockhis(age)=shocks+1
		    			end if
			    	else
				    	asum = asum + probsGoofy(sh)
					    sh=sh+1
    				end if
                end do
            else
    			wss=ran(iii)
	    		!	Assigning Markov productivity	!
		    	sstop=0
			    asum=0.0
    			s=1
	    		do while (sstop.eq.0)
		    		if (wss < asum + probTrue(history(age-1),s)) then
			    		history(age)=s
				    	sstop=1
					    if (s>states) then
    						print *,'Run #',j,' Person #',i,' Age =',age,' Markov =',s,' Random =',wss
	    					history(age)=states
		    			end if
    				else
	    				asum = asum + probTrue(history(age-1),s)
		    			s=s+1
			    	end if
    			end do

    			wss=ran(iii)
	    		!	Assigning transitory productivity	!
		    	sstop=0
			    asum=0.0
    			st=1
	    		do while (sstop.eq.0)
		    		if (wss < asum + probtTrue(st)) then
			    		transhis(age)=st
				    	sstop=1
    					if (st>trans) then
	    					print *,'Run #',j,' Person #',i,' Age =',age,' Transitory =',st,' Random =',wss
		    				transhis(age)=trans
			    		end if
				    else
					    asum = asum + probtTrue(st)
    					st=st+1
	    			end if
		    	end do

    			wss=ran(iii)
	    		!	Assigning expense shock	!
		    	sstop=0
			    asum=0.0
    			sh=1
	    		do while (sstop.eq.0)
		    		if (wss < asum + probsTrue(sh)) then
			    		shockhis(age)=sh
				    	sstop=1
					    if (sh>shocks+1) then
    						print *,'Run #',j,' Person #',i,' Age =',age,' Expense =',sh,' Random =',wss
	    					shockhis(age)=shocks+1
		    			end if
			    	else
				    	asum = asum + probsTrue(sh)
					    sh=sh+1
    				end if
                end do
            end if

            prior=(scorehis(age-1)-1.0)*scorebinsize
            posterior = (prior*probTrue(history(age-1),s)*probtTrue(st)*ProbsTrue(sh))/(prior*probTrue(history(age-1),s)*probtTrue(st)*ProbsTrue(sh) + (1-prior)*probGoofy(history(age-1),s)*probtGoofy(st)*ProbsGoofy(sh))
!            if (posterior>1.0) then
!                print *,'Posterior =', posterior
!                print *,'Prior =', prior
!                print *,'Transitory State is', st
!                pause
!            endif
            isprime = min(int(posterior/scorebinsize)+1,scorepoints)
            weight1 = (isprime*scorebinsize - posterior)/scorebinsize
            wss=ran(iii)
            if (wss<weight1) then
                scorehis(age) = isprime
            else
                scorehis(age) = min(isprime+1,scorepoints)
            end if

            if (goof) then
                avgscoreG(age)=avgscoreG(age)+scorehis(age)-1
            else
                avgscoreR(age)=avgscoreR(age)+scorehis(age)-1
            end if
            
			wage=phi(history(age))*A(transhis(age))*hbar(age)
			bin=bins(age,history(age),transhis(age))
			PopB(bin)=PopB(bin)+1
			outputB(bin)=outputB(bin)+wage

			ystig1=(stigma(age-1).eq.1)
			ystig2=(stigma(age-1).eq.2)
			ystig3=(stigma(age-1).eq.3)
			if (ystig1) then
				kpos(age)=bmax+1
			else if (ystig2) then
				kpos(age)=BRafter(age-1,history(age-1),transhis(age-1),scorehis(age-1),shockhis(age-1))
			else if (ystig3) then
				kpos(age)=DebtRolled(age-1,history(age-1),transhis(age-1),shockhis(age-1))
			else
				kpos(age)=BR(age-1,history(age-1),transhis(age-1),scorehis(age-1),kpos(age-1))
            end if
            preshock=kpos(age)
!			loss=nint(expense(history(age),shockhis(age))*gridfine)
			kpos(age)=max(kpos(age)-expenseloss(history(age),shockhis(age)),1)
			if (kpos(age) < bmax+1) then
				borrowers=borrowers+1.0
				avdoib=avdoib+((bmax+1.0-kpos(age))/gridfine)/wage
				borpostB(bin)=borpostB(bin)+1.0
				avdebincB(bin)=avdebincB(bin)+((bmax+1.0-kpos(age))/gridfine)/wage
				debtpostB(bin)=debtpostB(bin)+(bmax+1.0-kpos(age))/gridfine
				intinc=0.0
                if (goof) then
    				borrowersG=borrowersG+1.0
	    			avdoibG=avdoibG+((bmax+1.0-kpos(age))/gridfine)/wage
                else
    				borrowersR=borrowersR+1.0
	    			avdoibR=avdoibR+((bmax+1.0-kpos(age))/gridfine)/wage
                endif
			else
				intinc=((kpos(age)-bmax-1.0)/gridfine)*rs/(1+rs)
			end if
!			if (count<limitn) then
!				count=count+1
!				income(count)=wage+intinc
!				earnings(count)=wage
!				wealth(count)=(kpos(age)-bmax-1)/gridfine
!			endif

			if (wage>0.0) then
				output=output+wage
                if (goof) then
                    outputG=outputG+wage
                    incoG(age)=incoG(age)+wage
                else
                    outputR=outputR+wage
                    incoR(age)=incoR(age)+wage
                endif
				avlnincage(age)=avlnincage(age)+log(wage)
			else
				avlnincage(age)=avlnincage(age)-1.0d80
			end if

			run0=(V(age,history(age),transhis(age),scorehis(age),kpos(age))<VbarFS(age,history(age),transhis(age),scorehis(age)))
			defafter=(Vafter(age,history(age),transhis(age),scorehis(age),shockhis(age))<Vdef(age,history(age),transhis(age),scorehis(age),shockhis(age)))
			runedu=(Vedu(age,history(age),transhis(age),scorehis(age),kpos(age))<Vbaredu(age,history(age),transhis(age),scorehis(age)))
			defedu=(Vafteredu(age,history(age),transhis(age),scorehis(age),shockhis(age))<Vdefedu(age,history(age),transhis(age),scorehis(age),shockhis(age)))

			if (ystig1) then
				if (defafter) then
					defaults=defaults+1.0
					stigma(age)=3
					kprime=bmax+1		!	Rolled-over debts do not get counted
!					kprime=DebtRolled(age,history(age),transhis(age),shockhis(age))
					cons(age)=cons(age)+ConD(age,history(age),transhis(age),shockhis(age))+aut
                    welfarecum=welfarecum + betacum(age-1)*(Utility(aut+ConD(age,history(age),transhis(age),shockhis(age)),cigma)-coststig1)
					avlnconage(age)=avlnconage(age)+log(ConD(age,history(age),transhis(age),shockhis(age))+aut)
                    if (goof) then
    					avglnconG(age)=avglnconG(age)+log(ConD(age,history(age),transhis(age),shockhis(age))+aut)
                    else
    					avglnconR(age)=avglnconR(age)+log(ConD(age,history(age),transhis(age),shockhis(age))+aut)
                    end if
                    if (goof .and. not(defedu)) then
                        tooearly=tooearly+1.0
                        gain = gain + betacum(age-1)*(Vafteredu(age,history(age),transhis(age),scorehis(age),shockhis(age))-Vdefedu(age,history(age),transhis(age),scorehis(age),shockhis(age)))
                    end if
                    kprimeedu=bmax+1
				else
					stigma(age)=2
					kprime=BRafter(age,history(age),transhis(age),scorehis(age),shockhis(age))
					cons(age)=cons(age)+Con1(age,history(age),transhis(age),scorehis(age),shockhis(age))+aut
                    welfarecum=welfarecum + betacum(age-1)*(Utility(aut+Con1(age,history(age),transhis(age),scorehis(age),shockhis(age)),cigma)-coststig2)
					avlnconage(age)=avlnconage(age)+log(Con1(age,history(age),transhis(age),scorehis(age),shockhis(age))+aut)
                    if (goof) then
					    avglnconG(age)=avglnconG(age)+log(Con1(age,history(age),transhis(age),scorehis(age),shockhis(age))+aut)
                    else
					    avglnconR(age)=avglnconR(age)+log(Con1(age,history(age),transhis(age),scorehis(age),shockhis(age))+aut)
                    end if
                    if (goof) then
                        if (defedu) then
                            toolate=toolate+1.0
                            toolatebyage(age)=toolatebyage(age)+1.0
                            gain = gain - betacum(age-1)*(Vafteredu(age,history(age),transhis(age),scorehis(age),shockhis(age))-Vdefedu(age,history(age),transhis(age),scorehis(age),shockhis(age)))
                            overhang = overhang + q(age,history(age),scorehis(age),kprime)*max(bmax+1-kprime,0)/gridfine
                        else
        					kprimeedu=BRafteredu(age,history(age),transhis(age),scorehis(age),shockhis(age))
                            gain = gain + betacum(age-1)*(Vafteredu(age,history(age),transhis(age),scorehis(age),shockhis(age)) - (Vafter(age,history(age),transhis(age),scorehis(age),shockhis(age)) - beta(age)*Vprime(age,history(age),scorehis(age),kprime) + beta(age)*Vprimeedu(age,history(age),scorehis(age),kprime)))
                            if (kprime<min(kprimeedu,bmax+1)) then
                                overborrowing = overborrowing + (min(kprimeedu,bmax+1)-kprime)/gridfine
                                overborrq = overborrq + q(age,history(age),scorehis(age),kprime)*(min(kprimeedu,bmax+1)-kprime)/gridfine
                                countover = countover+1.0
                            elseif (kprimeedu<min(kprime,bmax+1)) then
                                countunder = countunder+1.0
                            end if
                        end if
                    end if
				end if
			elseif (run0) then
				bankruptcy=bankruptcy+1.0
                if (flag) then
                    recidivist=recidivist+1.0
                end if
                flag = .TRUE.
                if (wage>medianinc) then
                    richfiler = richfiler+1.0
                end if
				newbankruptcy(age)=newbankruptcy(age)+1.0
				filingsB(bin)=filingsB(bin)+1.0
				stigma(age)=1
				cons(age)=cons(age)+ConB(age,history(age),transhis(age))+aut
                welfarecum=welfarecum + betacum(age-1)*(Utility(aut+ConB(age,history(age),transhis(age)),cigma)-coststig1)
				avlnconage(age)=avlnconage(age)+log(ConB(age,history(age),transhis(age))+aut)
                if (goof) then
    				avglnconG(age)=avglnconG(age)+log(ConB(age,history(age),transhis(age))+aut)
                else
    				avglnconR(age)=avglnconR(age)+log(ConB(age,history(age),transhis(age))+aut)
                end if
				avdoid=avdoid+((bmax+1-kpos(age))/gridfine)/wage
				baddebt=baddebt+(bmax+1-kpos(age))/gridfine
				outdef=outdef+wage
                if (goof) then
    				bankruptcyG=bankruptcyG+1.0
    				newbankruptcyG(age)=newbankruptcyG(age)+1.0
    				avdoidG=avdoidG+((bmax+1-kpos(age))/gridfine)/wage
    				outdefG=outdefG+wage
                    if ((preshock<(bmax+1)) .and. not(ystig3)) then
                        profitsG=profitsG+(max((wage-expense(history(age),shockhis(age))-cbar)*garn,0.0))/(1+rb)
                        profitageG(age-1)=profitageG(age-1)+(max((wage-expense(history(age),shockhis(age))-cbar)*garn,0.0))/(1+rb)
                    end if
                    if (not(runedu)) then
                        gain = gain + betacum(age-1)*(Vedu(age,history(age),transhis(age),scorehis(age),kpos(age))-Vbaredu(age,history(age),transhis(age),scorehis(age)))
                        tooearly=tooearly+1.0
                    end if
                else
    				bankruptcyR=bankruptcyR+1.0
    				newbankruptcyR(age)=newbankruptcyR(age)+1.0
    				avdoidR=avdoidR+((bmax+1-kpos(age))/gridfine)/wage
    				outdefR=outdefR+wage
                    if ((preshock<(bmax+1)) .and. not(ystig3)) then
                        profitsR=profitsR+(max((wage-expense(history(age),shockhis(age))-cbar)*garn,0.0))/(1+rb)
                        profitageR(age-1)=profitageR(age-1)+(max((wage-expense(history(age),shockhis(age))-cbar)*garn,0.0))/(1+rb)
                    end if
                endif
				if (history(age)<history(age-1)) then
					increason=2		!	Persistent income shock	!
				else
					increason=1		!	No income shock	!
				endif
				if (transhis(age)<avtrans) then
					transreason=2		!	Transitory income shock	!
				else
					transreason=1		!	No income shock	!
				endif
				if (Vbaralt(age,history(age),transhis(age),scorehis(age))>VbarFS(age,history(age),transhis(age),scorehis(age))) then
					wannasave=wannasave+1.0
				endif

				BRbyreason(shockhis(age),increason,transreason)=BRbyreason(shockhis(age),increason,transreason)+1.0
                if (goof) then
    				BRbyreasonG(shockhis(age),increason,transreason)=BRbyreasonG(shockhis(age),increason,transreason)+1.0
                else
                    BRbyreasonR(shockhis(age),increason,transreason)=BRbyreasonR(shockhis(age),increason,transreason)+1.0
                end if
				kprime=bmax+1
                kprimeedu=bmax+1
			else
				stigma(age)=0
                if ((preshock<(bmax+1)) .and. not(ystig3)) then
                    if (goof) then
                        profitsG=profitsG+((bmax+1.0-preshock)/gridfine)/(1+rb)
                        profitageG(age-1)=profitageG(age-1)+((bmax+1.0-preshock)/gridfine)/(1+rb)
                    else
                        profitsR=profitsR+((bmax+1.0-preshock)/gridfine)/(1+rb)
                        profitageR(age-1)=profitageR(age-1)+((bmax+1.0-preshock)/gridfine)/(1+rb)
                    end if
                end if
				kprime=BR(age,history(age),transhis(age),scorehis(age),kpos(age))
				cons(age)=cons(age)+Con(age,history(age),transhis(age),scorehis(age),kpos(age))+aut
                welfarecum=welfarecum + betacum(age-1)*Utility(aut+Con(age,history(age),transhis(age),scorehis(age),kpos(age)),cigma)
				avlnconage(age)=avlnconage(age)+log(Con(age,history(age),transhis(age),scorehis(age),kpos(age))+aut)
                if (goof) then
    				avglnconG(age)=avglnconG(age)+log(Con(age,history(age),transhis(age),scorehis(age),kpos(age))+aut)
                else
    				avglnconR(age)=avglnconR(age)+log(Con(age,history(age),transhis(age),scorehis(age),kpos(age))+aut)
                end if
                if (goof) then
    				kprimeedu=BRedu(age,history(age),transhis(age),scorehis(age),kpos(age))
                    if (runedu) then
                        gain = gain - betacum(age-1)*(Vedu(age,history(age),transhis(age),scorehis(age),kpos(age))-Vbaredu(age,history(age),transhis(age),scorehis(age)))
                        toolate=toolate+1.0
                        toolatebyage(age)=toolatebyage(age)+1.0
                        overhang = overhang - (q(age,history(age),scorehis(age),kprime)*(kprime-1.0-bmax)/gridfine)
                    else
                        gain = gain + betacum(age-1)*(Vedu(age,history(age),transhis(age),scorehis(age),kpos(age)) - (V(age,history(age),transhis(age),scorehis(age),kpos(age)) - beta(age)*Vprime(age,history(age),scorehis(age),kprime) + beta(age)*Vprimeedu(age,history(age),scorehis(age),kprime)))
                        if (kprime<min(kprimeedu,bmax+1)) then
                            overborrowing = overborrowing + (min(kprimeedu,bmax+1)-kprime)/gridfine
                            overborrq = overborrq + q(age,history(age),scorehis(age),kprime)*(min(kprimeedu,bmax+1)-kprime)/gridfine
                            countover = countover+1.0
                        elseif (kprimeedu<min(kprime,bmax+1)) then
                            countunder = countunder+1.0
                        end if
                    end if
                end if
            end if

			if (kprime<bmax+1) then
				debtadd=-(q(age,history(age),scorehis(age),kprime)*(kprime-1.0-bmax)/gridfine)
				if (debtadd<0) then
					print *,'Debt<0'
					do kk=1,age
						print *,history(kk), transhis(kk), shockhis(kk)
					end do
					print *,'q=',q(age,history(age),scorehis(age),kprime)
					!pause
				else
					debt=debt+debtadd
					debtB(bin)=debtB(bin)+debtadd
					borrowB(bin)=borrowB(bin)+1.0
					debtbyage(age)=debtbyage(age)+debtadd
!					debtdist(age,history(age),transhis(age),shockhis(age))=debtdist(age,history(age),transhis(age),shockhis(age))+debtadd
					borrowante=borrowante+1.0
					borrow(age,history(age),transhis(age),shockhis(age))=borrow(age,history(age),transhis(age),shockhis(age)) + 1.0
					avborint=avborint-(1-q(age,history(age),scorehis(age),kprime))*(kprime-bmax-1)/gridfine
                    if (goof) then
        				debtG=debtG+debtadd
                        debtbyageG(age)=debtbyageG(age)+debtadd
    	    			borrowanteG=borrowanteG+1.0
    					borrowG(age,history(age),transhis(age),shockhis(age))=borrowG(age,history(age),transhis(age),shockhis(age)) + 1.0
    					avborintG=avborintG-(1-q(age,history(age),scorehis(age),kprime))*(kprime-bmax-1)/gridfine
    					interestageG(age)=interestageG(age)-(1-q(age,history(age),scorehis(age),kprime))*(kprime-bmax-1)/gridfine
                        profitsG=profitsG-debtadd
                        profitageG(age)=profitageG(age)-debtadd
!                        write (26,998) age,scorehis(age),debtadd,q(age,history(age),scorehis(age),kprime),q(age,history(age),1,kprime)
                    else
    				    debtR=debtR+debtadd
                        debtbyageR(age)=debtbyageR(age)+debtadd
        				borrowanteR=borrowanteR+1.0
    					borrowR(age,history(age),transhis(age),shockhis(age))=borrowR(age,history(age),transhis(age),shockhis(age)) + 1.0
    					avborintR=avborintR-(1-q(age,history(age),scorehis(age),kprime))*(kprime-bmax-1)/gridfine
    					interestageR(age)=interestageR(age)-(1-q(age,history(age),scorehis(age),kprime))*(kprime-bmax-1)/gridfine
                        profitsR=profitsR-debtadd
                        profitageR(age)=profitageR(age)-debtadd
!                        write (27,998) age,scorehis(age),debtadd,q(age,history(age),scorehis(age),kprime),q(age,history(age),scorepoints,kprime)
                    end if
				end if
			else
				debtadd=((kprime-1.0-bmax)/gridfine)/(1+rs)
				saving=saving+debtadd
				savebyage(age)=savebyage(age)+debtadd
			end if

		end do


		age=lifespan
		history(age)=history(age-1)
		transhis(age)=avtrans
		shockhis(age)=shocks+1
        scorehis(age)=scorehis(age-1)

		wage=phi(history(age))*hbar(age)
		ystig1=(stigma(age-1).eq.1)
		ystig2=(stigma(age-1).eq.2)
		ystig3=(stigma(age-1).eq.3)
		if (ystig1) then
			kpos(age)=bmax+1
            welfarecum=welfarecum + betacum(age-1)*(Utility(autlast+Con1(age,history(age),transhis(age),scorehis(age),shockhis(age)),cigma)-coststig2)
		else if (ystig2) then
			kpos(age)=BRafter(age-1,history(age-1),transhis(age-1),scorehis(age-1),shockhis(age-1))
            welfarecum=welfarecum + betacum(age-1)*Utility(autlast+Con(age,history(age),transhis(age),scorehis(age),kpos(age)),cigma)
		else if (ystig3) then
			kpos(age)=DebtRolled(age-1,history(age-1),transhis(age-1),shockhis(age-1))
            welfarecum=welfarecum + betacum(age-1)*(Utility(autlast+ConB(age,history(age),transhis(age)),cigma)-coststig1)
        else
			kpos(age)=BR(age-1,history(age-1),transhis(age-1),scorehis(age-1),kpos(age-1))
            welfarecum=welfarecum + betacum(age-1)*Utility(autlast+Con(age,history(age),transhis(age),scorehis(age),kpos(age)),cigma)
        end if
        if (kpos(age) < bmax+1) then
			borrowers=borrowers+1.0
			avdoib=avdoib+((bmax+1-kpos(age))/gridfine)/wage
            if (goof) then
        		borrowersG=borrowersG+1.0
    			avdoibG=avdoibG+((bmax+1-kpos(age))/gridfine)/wage
                if (not(ystig3)) then
                    profitsG=profitsG+((bmax+1-kpos(age))/gridfine)/(1+rb)
                    profitageG(age-1)=profitageG(age-1)+((bmax+1-kpos(age))/gridfine)/(1+rb)
                endif
            else
                borrowersR=borrowersR+1.0
            	avdoibR=avdoibR+((bmax+1-kpos(age))/gridfine)/wage
                if (not(ystig3)) then
                    profitsR=profitsR+((bmax+1-kpos(age))/gridfine)/(1+rb)
                    profitageR(age-1)=profitageR(age-1)+((bmax+1-kpos(age))/gridfine)/(1+rb)
                endif
            endif
			intinc=0.0
		else
			intinc=((kpos(age)-bmax-1)/gridfine)*rs/(1+rs)
		end if

		if (wage>0.0) then
			output=output+wage
            if (goof) then
                outputG=outputG+wage
                incoG(age)=incoG(age)+wage
            else
                outputR=outputR+wage
                incoR(age)=incoR(age)+wage
            endif
			avlnincage(age)=avlnincage(age)+log(wage)
		else
			avlnincage(age)=avlnincage(age)-1.0d80
		end if

		!	Note:	The only bankruptcies here are last periods' defaults		!
		if (ystig3) then
			bankruptcy=bankruptcy+1.0
            if (flag) then
                recidivist=recidivist+1.0
            end if
            if (wage>medianinc) then
                richfiler = richfiler+1.0
            end if
			newbankruptcy(age)=newbankruptcy(age)+1.0
			cons(age)=cons(age)+ConB(age,history(age),transhis(age))+autlast
			avlnconage(age)=avlnconage(age)+log(ConB(age,history(age),transhis(age))+autlast)
            if (goof) then
    			avglnconG(age)=avglnconG(age)+log(ConB(age,history(age),transhis(age))+autlast)
            else
    			avglnconR(age)=avglnconR(age)+log(ConB(age,history(age),transhis(age))+autlast)
            end if
			avdoid=avdoid+((bmax+1-kpos(age))/gridfine)/wage
			baddebt=baddebt+(bmax+1-kpos(age))/gridfine
			outdef=outdef+wage
            if (goof) then
        		bankruptcyG=bankruptcyG+1.0
				newbankruptcyG(age)=newbankruptcyG(age)+1.0
    			avdoidG=avdoidG+((bmax+1-kpos(age))/gridfine)/wage
				outdefG=outdefG+wage
            else
        		bankruptcyR=bankruptcyR+1.0
				newbankruptcyR(age)=newbankruptcyR(age)+1.0
    			avdoidR=avdoidR+((bmax+1-kpos(age))/gridfine)/wage
				outdefR=outdefR+wage
            endif
			increason=1			!	No income shock	!
			transreason=1		!	No income shock	!
			BRbyreason(shockhis(age),increason,transreason)=BRbyreason(shockhis(age),increason,transreason)+1.0
            if (goof) then
                BRbyreasonG(shockhis(age),increason,transreason)=BRbyreasonG(shockhis(age),increason,transreason)+1.0
            else
                BRbyreasonR(shockhis(age),increason,transreason)=BRbyreasonR(shockhis(age),increason,transreason)+1.0
            end if
		else if (ystig1) then
			cons(age)=cons(age)+Con1(age,history(age),transhis(age),scorehis(age),shockhis(age))+autlast
			avlnconage(age)=avlnconage(age)+log(Con1(age,history(age),transhis(age),scorehis(age),shockhis(age))+autlast)
            if (goof) then
    			avglnconG(age)=avglnconG(age)+log(Con1(age,history(age),transhis(age),scorehis(age),shockhis(age))+autlast)
            else
    			avglnconR(age)=avglnconR(age)+log(Con1(age,history(age),transhis(age),scorehis(age),shockhis(age))+autlast)
            end if
		else
			cons(age)=cons(age)+Con(age,history(age),transhis(age),scorehis(age),kpos(age))+autlast
			avlnconage(age)=avlnconage(age)+log(Con(age,history(age),transhis(age),scorehis(age),kpos(age))+autlast)
            if (goof) then
    			avglnconG(age)=avglnconG(age)+log(Con(age,history(age),transhis(age),scorehis(age),kpos(age))+autlast)
            else
    			avglnconR(age)=avglnconR(age)+log(Con(age,history(age),transhis(age),scorehis(age),kpos(age))+autlast)
            end if
        end if
        realizedEV = realizedEV + welfarecum
        if (goof) then
            realizedG=realizedG+welfarecum
        else
            realizedR=realizedR+welfarecum
        end if
        if (welfarecum>bestluck) then
            bestluck = welfarecum
        end if
        if (welfarecum<worstluck) then
            worstluck = welfarecum
        end if
        
	end do

    realizedwelfare(j) = realizedEV/Pop
    bankrtotal(j)=(bankruptcy/Pop)/lifespan
    netprofitG(j)=(profitsG/Pop)/lifespan
    netprofitR(j)=(profitsR/Pop)/lifespan
    if (countR>0.0) then
        welfareR(j)=realizedR/countR
        bankrR(j)=(bankruptcyR/countR)/lifespan
    end if
    if (countG>0.0) then
        welfareG(j)=realizedG/countG
        bankrG(j)=(bankruptcyG/countG)/lifespan
        overborrowtotal(j)=(overborrowing/countG)/lifespan
        overborrqtotal(j)=(overborrq/countG)/lifespan
        overhangtotal(j)=(overhang/countG)/lifespan
        gaintotal(j)=(gain/countG)/lifespan
        totalover(j)=(countover/countG)/(lifespan-1)
        totalunder(j)=(countunder/countG)/(lifespan-1)
        totalearly(j)=(tooearly/countG)/(lifespan-1)
        totallate(j)=(toolate/countG)/(lifespan-1)
    end if
	defaulttotal(j)=(defaults/Pop)/lifespan
	if (debt>0.0) then
		avinttotal(j)=avborint/debt
		debttotal(j)=(debt/Pop)/lifespan
		borrowerantetotal(j)=(borrowante/Pop)/lifespan
	else
		avinttotal(j)=0.0
		debttotal(j)=0.0
		borrowerantetotal(j)=0.0
	end if
	if (debtG>0.0) then
		avintG(j)=avborintG/debtG
		debttotalG(j)=(debtG/countG)/lifespan
		borroweranteG(j)=(borrowanteG/countG)/lifespan
	else
		avintG(j)=0.0
		debttotalG(j)=0.0
		borroweranteG(j)=0.0
	end if
	if (debtR>0.0) then
		avintR(j)=avborintR/debtR
		debttotalR(j)=(debtR/countR)/lifespan
		borroweranteR(j)=(borrowanteR/countR)/lifespan
	else
		avintR(j)=0.0
		debttotalR(j)=0.0
		borroweranteR(j)=0.0
	end if
	savetotal(j)=(saving/Pop)/lifespan
	if (borrowers>0.0) then
		borrowerposttotal(j)=(borrowers/Pop)/lifespan
		debtincbortotal(j)=avdoib/borrowers
	else
		borrowerposttotal(j)=0.0
		debtincbortotal(j)=0.0
	end if
	if (borrowersG>0.0) then
		borrowerpostG(j)=(borrowersG/countG)/lifespan
		debtincborG(j)=avdoibG/borrowersG
	else
		borrowerpostG(j)=0.0
		debtincborG(j)=0.0
	end if
	if (borrowersR>0.0) then
		borrowerpostR(j)=(borrowersR/countR)/lifespan
		debtincborR(j)=avdoibR/borrowersR
	else
		borrowerpostR(j)=0.0
		debtincborR(j)=0.0
	end if
	if (bankruptcy>0.0) then
		debtincdeftotal(j)=avdoid/bankruptcy
		writeoff(j)=baddebt/bankruptcy
		savers(j)=wannasave/bankruptcy
        repeatfrac(j)=recidivist/bankruptcy
        richfrac(j)=richfiler/bankruptcy
	else
		debtincdeftotal(j)=0.0
		writeoff(j)=0.0
		savers(j)=0.0
        repeatfrac(j)=0.0
        richfrac(j)=0.0
	end if
	if (bankruptcyG>0.0) then
		debtincdefG(j)=avdoidG/bankruptcyG
	else
		debtincdefG(j)=0.0
	end if
	if (bankruptcyR>0.0) then
		debtincdefR(j)=avdoidR/bankruptcyR
	else
		debtincdefR(j)=0.0
	end if

end do

avlnconage = (avlnconage/Pop)/Monte1
avglnconG = (avglnconG/Pop)/Monte1
avglnconR = (avglnconR/Pop)/Monte1
avlnincage = (avlnincage/Pop)/Monte1
avlncon=0.0
avlninc=0.0
do age=1,lifespan
	avlncon = avlncon + avlnconage(age)
	avlninc = avlninc + avlnincage(age)
end do
avlncon = avlncon/lifespan
avlninc = avlninc/lifespan
avborint=0.0
debt=0.0
do i=1,Monte1
    avborint=avborint+avinttotal(i)*debttotal(i)
    debt = debt + debttotal(i)
end do
avborint1=avborint/debt		!	Debt-weighted across Monte1 simulations!

varlnconage = 0.0
varlnconG = 0.0
varlnconR = 0.0
varlnincage = 0.0
varint=0.0

do j=(Monte1+1),Monte
	debt=0.0
	debtG=0.0
	debtR=0.0
	saving=0.0
	avborint=0.0
	avborintG=0.0
	avborintR=0.0
	borrowante=0.0
	borrowanteG=0.0
	borrowanteR=0.0
	bankruptcy=0.0
    recidivist=0.0
    richfiler=0.0
	bankruptcyG=0.0
	bankruptcyR=0.0
	wannasave=0.0
	defaults=0.0
	borrowers=0.0
	borrowersG=0.0
	borrowersR=0.0
	avdoib=0.0
	avdoibG=0.0
	avdoibR=0.0
	avdoid=0.0
	avdoidG=0.0
	avdoidR=0.0
	baddebt=0.0
    realizedEV=0.0
    realizedG=0.0
    realizedR=0.0
    overborrowing=0.0
    overborrq=0.0
    overhang=0.0
    tooearly=0.0
    toolate=0.0
    gain=0.0
    profitsG=0.0
    profitsR=0.0
    countG = 0.0
    countR = 0.0
    countover = 0.0
    countunder = 0.0
    print *,'Monte=',j

	do i=1,Pop
        welfarecum=0.0
        flag=.FALSE.
        age=1
		wss=ran(iii)
		!	Assigning type	!
        goof=(wss<lambda)
        if (goof) then
            PopG = PopG + 1
            countG=countG+1.0
        else
            PopR = PopR + 1
            countR=countR+1.0
        end if

		wss=ran(iii)
		!	Assigning Markov productivity	!
		sstop=0
		asum=0.0
		s=1
		do while (sstop.eq.0)
			if (wss < asum + prob0(s)) then
				history(1)=s
				sstop=1
				if (s>states) then
					print *,'Run #',j,' Markov =',s,' Random =',wss
					history(1)=states
				end if
			else
				asum = asum + prob0(s)
				s=s+1
			end if
		end do

		wss=ran(iii)
		!	Assigning transitory productivity	!
		sstop=0
		asum=0.0
		st=1
        if (goof) then
    		do while (sstop.eq.0)
	    		if (wss < asum + probtGoofy(st)) then
    				transhis(1)=st
	    			sstop=1
		    		if (st>trans) then
			    		print *,'Run #',j,' Transitory =',st,' Random =',wss
				    	transhis(1)=trans
    				end if
	    		else
		    		asum = asum + probtGoofy(st)
    				st=st+1
	    		end if
            end do
        else   
		    do while (sstop.eq.0)
			    if (wss < asum + probtTrue(st)) then
				    transhis(1)=st
    				sstop=1
	    			if (st>trans) then
		    			print *,'Run #',j,' Transitory =',st,' Random =',wss
			    		transhis(1)=trans
				    end if
    			else
	    			asum = asum + probtTrue(st)
		    		st=st+1
    			end if
            end do
        end if
        posterior = ((1-lambda)*probtTrue(st))/((1-lambda)*probtTrue(st) + lambda*probtGoofy(st))
        isprime = min(int(posterior/scorebinsize)+1,scorepoints)
        weight1 = (isprime*scorebinsize - posterior)/scorebinsize
        wss=ran(iii)
        if (wss<weight1) then
            scorehis(1) = isprime
        else
            scorehis(1) = min(isprime+1,scorepoints)
        end if
        
        if (goof) then
            avgscoreG(1)=avgscoreG(1)+scorehis(1)-1
        else
            avgscoreR(1)=avgscoreR(1)+scorehis(1)-1
        end if

		shockhis(1)=shocks+1
		kpos(1)=bmax+1

		wage=phi(history(1))*A(transhis(1))*hbar(1)
		bin=bins(1,history(1),transhis(1))
		PopB(bin)=PopB(bin)+1
		outputB(bin)=outputB(bin)+wage

		if (wage>0.0) then
			output=output+wage
            if (goof) then
                outputG=outputG+wage
                incoG(1)=incoG(1)+wage
            else
                outputR=outputR+wage
                incoR(1)=incoR(1)+wage
            endif
			varlnincage(1)=varlnincage(1)+(log(wage)-avlnincage(1))**2
		else
			varlnincage(1)=varlnincage(1)+(-1.0d80-avlnincage(1))**2
		end if

        stigma(1)=0
    	kprime=BR(1,history(1),transhis(1),scorehis(1),kpos(1))
		cons(1)=cons(1)+Con(1,history(1),transhis(1),scorehis(1),kpos(1))+aut
        welfarecum=welfarecum + Utility(aut+Con(1,history(1),transhis(1),scorehis(1),kpos(1)),cigma)
    	varlnconage(1)=varlnconage(1)+(log(Con(1,history(1),transhis(1),scorehis(1),kpos(1))+aut)-avlnconage(1))**2
        if (goof) then
        	varlnconG(1)=varlnconG(1)+(log(Con(1,history(1),transhis(1),scorehis(1),kpos(1))+aut)-avglnconG(1))**2
        else
        	varlnconR(1)=varlnconR(1)+(log(Con(1,history(1),transhis(1),scorehis(1),kpos(1))+aut)-avglnconR(1))**2
        end if

        if (goof) then
            kprimeedu=BRedu(1,history(1),transhis(1),scorehis(1),kpos(1))
            gain = gain + Vedu(1,history(1),transhis(1),scorehis(1),kpos(1)) - (V(1,history(1),transhis(1),scorehis(1),kpos(1)) - beta(1)*Vprime(1,history(1),scorehis(1),kprime) + beta(1)*Vprimeedu(1,history(1),scorehis(1),kprime))
            if (kprime<min(kprimeedu,bmax+1)) then
                overborrowing = overborrowing + (min(kprimeedu,bmax+1)-kprime)/gridfine
                overborrq = overborrq + q(1,history(1),scorehis(1),kprime)*(min(kprimeedu,bmax+1)-kprime)/gridfine
                countover=countover+1.0
            elseif (kprimeedu<min(kprime,bmax+1)) then
                countunder=countunder+1.0
            end if
        end if

		if (kprime<bmax+1) then
			debtadd=-(q(1,history(1),scorehis(1),kprime)*(kprime-1-bmax)/gridfine)
			if (debtadd<0) then
				print *,'Debt<0'
				print *,'q=',q(1,history(1),scorehis(1),kprime)
				!pause
			else
				debt=debt+debtadd
				debtB(bin)=debtB(bin)+debtadd
				borrowB(bin)=borrowB(bin)+1.0
				debtbyage(1)=debtbyage(1)+debtadd
				borrowante=borrowante+1.0
				borrow(1,history(1),transhis(1),shockhis(1))=borrow(1,history(1),transhis(1),shockhis(1)) + 1.0
				avborint=avborint-(1-q(1,history(1),scorehis(1),kprime))*(kprime-bmax-1)/gridfine
                varint = varint + debtadd*(1/q(1,history(1),scorehis(1),kprime) -1 -avborint1)**2
                if (goof) then
    				debtG=debtG+debtadd
                    debtbyageG(1)=debtbyageG(1)+debtadd
    				borrowanteG=borrowanteG+1.0
    				borrowG(1,history(1),transhis(1),shockhis(1))=borrowG(1,history(1),transhis(1),shockhis(1)) + 1.0
    				avborintG=avborintG-(1-q(1,history(1),scorehis(1),kprime))*(kprime-bmax-1)/gridfine
    				interestageG(1)=interestageG(1)-(1-q(1,history(1),scorehis(1),kprime))*(kprime-bmax-1)/gridfine
                    profitsG=profitsG-debtadd
                    profitageG(1)=profitageG(1)-debtadd
!                    write (26,998) age,scorehis(1),debtadd,q(1,history(1),scorehis(1),kprime),q(1,history(1),1,kprime)
                else
    				debtR=debtR+debtadd
                    debtbyageR(1)=debtbyageR(1)+debtadd
    				borrowanteR=borrowanteR+1.0
    				borrowR(1,history(1),transhis(1),shockhis(1))=borrowR(1,history(1),transhis(1),shockhis(1)) + 1.0
    				avborintR=avborintR-(1-q(1,history(1),scorehis(1),kprime))*(kprime-bmax-1)/gridfine
    				interestageR(1)=interestageR(1)-(1-q(1,history(1),scorehis(1),kprime))*(kprime-bmax-1)/gridfine
                    profitsR=profitsR-debtadd
                    profitageR(1)=profitageR(1)-debtadd
!                    write (27,998) age,scorehis(1),debtadd,q(1,history(1),scorehis(1),kprime),q(1,history(1),scorepoints,kprime)
                end if
			end if
		else
			debtadd=((kprime-1-bmax)/gridfine)/(1+rs)
			saving=saving+debtadd
			savebyage(1)=savebyage(1)+debtadd
		end if


		do age=2,lifespan-1
            if (goof) then
    			wss=ran(iii)
	    		!	Assigning Markov productivity	!
		    	sstop=0
			    asum=0.0
    			s=1
	    		do while (sstop.eq.0)
		    		if (wss < asum + probGoofy(history(age-1),s)) then
			    		history(age)=s
				    	sstop=1
					    if (s>states) then
    						print *,'Run #',j,' Person #',i,' Age =',age,' Markov =',s,' Random =',wss
	    					history(age)=states
		    			end if
    				else
	    				asum = asum + probGoofy(history(age-1),s)
		    			s=s+1
			    	end if
    			end do

    			wss=ran(iii)
	    		!	Assigning transitory productivity	!
		    	sstop=0
			    asum=0.0
    			st=1
	    		do while (sstop.eq.0)
		    		if (wss < asum + probtGoofy(st)) then
			    		transhis(age)=st
				    	sstop=1
    					if (st>trans) then
	    					print *,'Run #',j,' Person #',i,' Age =',age,' Transitory =',st,' Random =',wss
		    				transhis(age)=trans
			    		end if
				    else
					    asum = asum + probtGoofy(st)
    					st=st+1
	    			end if
		    	end do

    			wss=ran(iii)
	    		!	Assigning expense shock	!
		    	sstop=0
			    asum=0.0
    			sh=1
	    		do while (sstop.eq.0)
		    		if (wss < asum + probsGoofy(sh)) then
			    		shockhis(age)=sh
				    	sstop=1
					    if (sh>shocks+1) then
    						print *,'Run #',j,' Person #',i,' Age =',age,' Expense =',sh,' Random =',wss
	    					shockhis(age)=shocks+1
		    			end if
			    	else
				    	asum = asum + probsGoofy(sh)
					    sh=sh+1
    				end if
                end do
                
            else
    			wss=ran(iii)
	    		!	Assigning Markov productivity	!
		    	sstop=0
			    asum=0.0
    			s=1
	    		do while (sstop.eq.0)
		    		if (wss < asum + probTrue(history(age-1),s)) then
			    		history(age)=s
				    	sstop=1
					    if (s>states) then
    						print *,'Run #',j,' Person #',i,' Age =',age,' Markov =',s,' Random =',wss
	    					history(age)=states
		    			end if
    				else
	    				asum = asum + probTrue(history(age-1),s)
		    			s=s+1
			    	end if
    			end do

    			wss=ran(iii)
	    		!	Assigning transitory productivity	!
		    	sstop=0
			    asum=0.0
    			st=1
	    		do while (sstop.eq.0)
		    		if (wss < asum + probtTrue(st)) then
			    		transhis(age)=st
				    	sstop=1
    					if (st>trans) then
	    					print *,'Run #',j,' Person #',i,' Age =',age,' Transitory =',st,' Random =',wss
		    				transhis(age)=trans
			    		end if
				    else
					    asum = asum + probtTrue(st)
    					st=st+1
	    			end if
		    	end do

    			wss=ran(iii)
	    		!	Assigning expense shock	!
		    	sstop=0
			    asum=0.0
    			sh=1
	    		do while (sstop.eq.0)
		    		if (wss < asum + probsTrue(sh)) then
			    		shockhis(age)=sh
				    	sstop=1
					    if (sh>shocks+1) then
    						print *,'Run #',j,' Person #',i,' Age =',age,' Expense =',sh,' Random =',wss
	    					shockhis(age)=shocks+1
		    			end if
			    	else
				    	asum = asum + probsTrue(sh)
					    sh=sh+1
    				end if
                end do
            end if

            prior=(scorehis(age-1)-1.0)*scorebinsize
            posterior = (prior*probTrue(history(age-1),s)*probtTrue(st)*ProbsTrue(sh))/(prior*probTrue(history(age-1),s)*probtTrue(st)*ProbsTrue(sh) + (1-prior)*probGoofy(history(age-1),s)*probtGoofy(st)*ProbsGoofy(sh))
            isprime = min(int(posterior/scorebinsize)+1,scorepoints)
            weight1 = (isprime*scorebinsize - posterior)/scorebinsize
            wss=ran(iii)
            if (wss<weight1) then
                scorehis(age) = isprime
            else
                scorehis(age) = min(isprime+1,scorepoints)
            end if

            if (goof) then
                avgscoreG(age)=avgscoreG(age)+scorehis(age)-1
            else
                avgscoreR(age)=avgscoreR(age)+scorehis(age)-1
            end if

			wage=phi(history(age))*A(transhis(age))*hbar(age)
			bin=bins(age,history(age),transhis(age))
			PopB(bin)=PopB(bin)+1
			outputB(bin)=outputB(bin)+wage

			ystig1=(stigma(age-1).eq.1)
			ystig2=(stigma(age-1).eq.2)
			ystig3=(stigma(age-1).eq.3)
			if (ystig1) then
				kpos(age)=bmax+1
			else if (ystig2) then
				kpos(age)=BRafter(age-1,history(age-1),transhis(age-1),scorehis(age-1),shockhis(age-1))
			else if (ystig3) then
				kpos(age)=DebtRolled(age-1,history(age-1),transhis(age-1),shockhis(age-1))
			else
				kpos(age)=BR(age-1,history(age-1),transhis(age-1),scorehis(age-1),kpos(age-1))
			end if
            preshock=kpos(age)
!			loss=nint(expense(history(age),shockhis(age))*gridfine)
			kpos(age)=max(kpos(age)-expenseloss(history(age),shockhis(age)),1)
			if (kpos(age) < bmax+1) then
				borrowers=borrowers+1.0
				avdoib=avdoib+((bmax+1-kpos(age))/gridfine)/wage
				borpostB(bin)=borpostB(bin)+1.0
				avdebincB(bin)=avdebincB(bin)+((bmax+1-kpos(age))/gridfine)/wage
				debtpostB(bin)=debtpostB(bin)+(bmax+1-kpos(age))/gridfine
                if (goof) then
    				borrowersG=borrowersG+1.0
	    			avdoibG=avdoibG+((bmax+1-kpos(age))/gridfine)/wage
                else
    				borrowersR=borrowersR+1.0
	    			avdoibR=avdoibR+((bmax+1-kpos(age))/gridfine)/wage
                endif
			end if

			if (wage>0.0) then
				output=output+wage
                if (goof) then
                    outputG=outputG+wage
                    incoG(age)=incoG(age)+wage
                else
                    outputR=outputR+wage
                    incoR(age)=incoR(age)+wage
                endif
				varlnincage(age)=varlnincage(age)+(log(wage)-avlnincage(age))**2
			else
				varlnincage(age)=varlnincage(age)+(-1.0d80-avlnincage(age))**2
			end if

			run0=(V(age,history(age),transhis(age),scorehis(age),kpos(age))<VbarFS(age,history(age),transhis(age),scorehis(age)))
			defafter=(Vafter(age,history(age),transhis(age),scorehis(age),shockhis(age))<Vdef(age,history(age),transhis(age),scorehis(age),shockhis(age)))
			runedu=(Vedu(age,history(age),transhis(age),scorehis(age),kpos(age))<Vbaredu(age,history(age),transhis(age),scorehis(age)))
			defedu=(Vafteredu(age,history(age),transhis(age),scorehis(age),shockhis(age))<Vdefedu(age,history(age),transhis(age),scorehis(age),shockhis(age)))

			if (ystig1) then
				if (defafter) then
					defaults=defaults+1.0
					stigma(age)=3
					kprime=bmax+1		!	Rolled-over debts do not get counted
!					kprime=DebtRolled(age,history(age),transhis(age),shockhis(age))
					cons(age)=cons(age)+ConD(age,history(age),transhis(age),shockhis(age))+aut
                    welfarecum=welfarecum + betacum(age-1)*(Utility(aut+ConD(age,history(age),transhis(age),shockhis(age)),cigma)-coststig1)
					varlnconage(age)=varlnconage(age)+(log(ConD(age,history(age),transhis(age),shockhis(age))+aut)-avlnconage(age))**2
                    if (goof) then
    					varlnconG(age)=varlnconG(age)+(log(ConD(age,history(age),transhis(age),shockhis(age))+aut)-avglnconG(age))**2
                    else
    					varlnconR(age)=varlnconR(age)+(log(ConD(age,history(age),transhis(age),shockhis(age))+aut)-avglnconR(age))**2
                    end if
                    if (goof .and. not(defedu)) then
                        tooearly=tooearly+1.0
                        gain = gain + betacum(age-1)*(Vafteredu(age,history(age),transhis(age),scorehis(age),shockhis(age))-Vdefedu(age,history(age),transhis(age),scorehis(age),shockhis(age)))
                    end if
                    kprimeedu=bmax+1
				else
					stigma(age)=2
					kprime=BRafter(age,history(age),transhis(age),scorehis(age),shockhis(age))
					cons(age)=cons(age)+Con1(age,history(age),transhis(age),scorehis(age),shockhis(age))+aut
                    welfarecum=welfarecum + betacum(age-1)*(Utility(aut+Con1(age,history(age),transhis(age),scorehis(age),shockhis(age)),cigma)-coststig2)
					varlnconage(age)=varlnconage(age)+(log(Con1(age,history(age),transhis(age),scorehis(age),shockhis(age))+aut)-avlnconage(age))**2
                    if (goof) then
    					varlnconG(age)=varlnconG(age)+(log(Con1(age,history(age),transhis(age),scorehis(age),shockhis(age))+aut)-avglnconG(age))**2
                    else
    					varlnconR(age)=varlnconR(age)+(log(Con1(age,history(age),transhis(age),scorehis(age),shockhis(age))+aut)-avglnconR(age))**2
                    end if
                    if (goof) then
                        if (defedu) then
                            toolate=toolate+1.0
                            toolatebyage(age)=toolatebyage(age)+1.0
                            gain = gain - betacum(age-1)*(Vafteredu(age,history(age),transhis(age),scorehis(age),shockhis(age))-Vdefedu(age,history(age),transhis(age),scorehis(age),shockhis(age)))
                            overhang = overhang + q(age,history(age),scorehis(age),kprime)*max(bmax+1-kprime,0)/gridfine
                        else
        					kprimeedu=BRafteredu(age,history(age),transhis(age),scorehis(age),shockhis(age))
                            gain = gain + betacum(age-1)*(Vafteredu(age,history(age),transhis(age),scorehis(age),shockhis(age)) - (Vafter(age,history(age),transhis(age),scorehis(age),shockhis(age)) - beta(age)*Vprime(age,history(age),scorehis(age),kprime) + beta(age)*Vprimeedu(age,history(age),scorehis(age),kprime)))
                            if (kprime<min(kprimeedu,bmax+1)) then
                                overborrowing = overborrowing + (min(kprimeedu,bmax+1)-kprime)/gridfine
                                overborrq = overborrq + q(age,history(age),scorehis(age),kprime)*(min(kprimeedu,bmax+1)-kprime)/gridfine
                                countover = countover+1.0
                            elseif (kprimeedu<min(kprime,bmax+1)) then
                                countunder = countunder+1.0
                            end if
                        end if
                    end if
				end if
			else if (run0) then
				bankruptcy=bankruptcy+1.0
                if (flag) then
                    recidivist=recidivist+1.0
                end if
                flag = .TRUE.
                if (wage>medianinc) then
                    richfiler = richfiler+1.0
                end if
				newbankruptcy(age)=newbankruptcy(age)+1.0
				filingsB(bin)=filingsB(bin)+1.0
				stigma(age)=1
				cons(age)=cons(age)+ConB(age,history(age),transhis(age))+aut
                welfarecum=welfarecum + betacum(age-1)*(Utility(aut+ConB(age,history(age),transhis(age)),cigma)-coststig1)
				varlnconage(age)=varlnconage(age)+(log(ConB(age,history(age),transhis(age))+aut)-avlnconage(age))**2
				avdoid=avdoid+((bmax+1-kpos(age))/gridfine)/wage
				baddebt=baddebt+(bmax+1-kpos(age))/gridfine
				outdef=outdef+wage
                if (goof) then
    				varlnconG(age)=varlnconG(age)+(log(ConB(age,history(age),transhis(age))+aut)-avglnconG(age))**2
    				bankruptcyG=bankruptcyG+1.0
    				newbankruptcyG(age)=newbankruptcyG(age)+1.0
    				avdoidG=avdoidG+((bmax+1-kpos(age))/gridfine)/wage
    				outdefG=outdefG+wage
                    if ((preshock<(bmax+1)) .and. not(ystig3)) then
                        profitsG=profitsG+(max((wage-expense(history(age),shockhis(age))-cbar)*garn,0.0))/(1+rb)
                        profitageG(age-1)=profitageG(age-1)+(max((wage-expense(history(age),shockhis(age))-cbar)*garn,0.0))/(1+rb)
                    end if
                    if (not(runedu)) then
                        tooearly=tooearly+1.0
                        gain = gain + betacum(age-1)*(Vedu(age,history(age),transhis(age),scorehis(age),kpos(age))-Vbaredu(age,history(age),transhis(age),scorehis(age)))
                    end if
                else
    				varlnconR(age)=varlnconR(age)+(log(ConB(age,history(age),transhis(age))+aut)-avglnconR(age))**2
    				bankruptcyR=bankruptcyR+1.0
    				newbankruptcyR(age)=newbankruptcyR(age)+1.0
    				avdoidR=avdoidR+((bmax+1-kpos(age))/gridfine)/wage
    				outdefR=outdefR+wage
                    if ((preshock<(bmax+1)) .and. not(ystig3)) then
                        profitsR=profitsR+(max((wage-expense(history(age),shockhis(age))-cbar)*garn,0.0))/(1+rb)
                        profitageR(age-1)=profitageR(age-1)+(max((wage-expense(history(age),shockhis(age))-cbar)*garn,0.0))/(1+rb)
                    end if
                endif
				if (history(age)<history(age-1)) then
					increason=2		!	Persistent income shock	!
				else
					increason=1		!	No income shock	!
				endif
				if (transhis(age)<avtrans) then
					transreason=2		!	Transitory income shock	!
				else
					transreason=1		!	No income shock	!
				endif
				if (Vbaralt(age,history(age),transhis(age),scorehis(age))>VbarFS(age,history(age),transhis(age),scorehis(age))) then
					wannasave=wannasave+1.0
				endif

				BRbyreason(shockhis(age),increason,transreason)=BRbyreason(shockhis(age),increason,transreason)+1.0
                if (goof) then
    				BRbyreasonG(shockhis(age),increason,transreason)=BRbyreasonG(shockhis(age),increason,transreason)+1.0
                else
                    BRbyreasonR(shockhis(age),increason,transreason)=BRbyreasonR(shockhis(age),increason,transreason)+1.0
                end if
				kprime=bmax+1
                kprimeedu=bmax+1
			else
				stigma(age)=0
                if ((preshock<(bmax+1)) .and. not(ystig3)) then
                    if (goof) then
                        profitsG=profitsG+((bmax+1.0-preshock)/gridfine)/(1+rb)
                        profitageG(age-1)=profitageG(age-1)+((bmax+1.0-preshock)/gridfine)/(1+rb)
                    else
                        profitsR=profitsR+((bmax+1.0-preshock)/gridfine)/(1+rb)
                        profitageR(age-1)=profitageR(age-1)+((bmax+1.0-preshock)/gridfine)/(1+rb)
                    end if
                end if
				kprime=BR(age,history(age),transhis(age),scorehis(age),kpos(age))
				cons(age)=cons(age)+Con(age,history(age),transhis(age),scorehis(age),kpos(age))+aut
                welfarecum=welfarecum + betacum(age-1)*Utility(aut+Con(age,history(age),transhis(age),scorehis(age),kpos(age)),cigma)
				varlnconage(age)=varlnconage(age)+(log(Con(age,history(age),transhis(age),scorehis(age),kpos(age))+aut)-avlnconage(age))**2
                if (goof) then
    				varlnconG(age)=varlnconG(age)+(log(Con(age,history(age),transhis(age),scorehis(age),kpos(age))+aut)-avglnconG(age))**2
                else
    				varlnconR(age)=varlnconR(age)+(log(Con(age,history(age),transhis(age),scorehis(age),kpos(age))+aut)-avglnconR(age))**2
                end if
                if (goof) then
    				kprimeedu=BRedu(age,history(age),transhis(age),scorehis(age),kpos(age))
                    if (runedu) then
                        gain = gain - betacum(age-1)*(Vedu(age,history(age),transhis(age),scorehis(age),kpos(age))-Vbaredu(age,history(age),transhis(age),scorehis(age)))
                        toolate=toolate+1.0
                        toolatebyage(age)=toolatebyage(age)+1.0
                        overhang = overhang - (q(age,history(age),scorehis(age),kprime)*(kprime-1.0-bmax)/gridfine)
                    else
                        gain = gain + betacum(age-1)*(Vedu(age,history(age),transhis(age),scorehis(age),kpos(age)) - (V(age,history(age),transhis(age),scorehis(age),kpos(age)) - beta(age)*Vprime(age,history(age),scorehis(age),kprime) + beta(age)*Vprimeedu(age,history(age),scorehis(age),kprime)))
                        if (kprime<min(kprimeedu,bmax+1)) then
                            overborrowing = overborrowing + (min(kprimeedu,bmax+1)-kprime)/gridfine
                            overborrq = overborrq + q(age,history(age),scorehis(age),kprime)*(min(kprimeedu,bmax+1)-kprime)/gridfine
                            countover = countover+1.0
                        elseif (kprimeedu<min(kprime,bmax+1)) then
                            countunder = countunder+1.0
                        end if
                    end if
                end if
			end if

			if (kprime<bmax+1) then
				debtadd=-(q(age,history(age),scorehis(age),kprime)*(kprime-1.0-bmax)/gridfine)
				if (debtadd<0) then
					print *,'Debt<0'
					do kk=1,age
						print *,history(kk), transhis(kk), shockhis(kk)
					end do
					print *,'q=',q(age,history(age),scorehis(age),kprime)
					!pause
				else
					debt=debt+debtadd
					debtB(bin)=debtB(bin)+debtadd
					borrowB(bin)=borrowB(bin)+1.0
					debtbyage(age)=debtbyage(age)+debtadd
!					debtdist(age,history(age),transhis(age),shockhis(age))=debtdist(age,history(age),transhis(age),shockhis(age))+debtadd
					borrowante=borrowante+1.0
					borrow(age,history(age),transhis(age),shockhis(age))=borrow(age,history(age),transhis(age),shockhis(age)) + 1.0
					avborint=avborint-(1-q(age,history(age),scorehis(age),kprime))*(kprime-bmax-1)/gridfine
                    varint = varint + debtadd*(1/q(age,history(age),scorehis(age),kprime) -1 -avborint1)**2
                    if (goof) then
        				debtG=debtG+debtadd
                        debtbyageG(age)=debtbyageG(age)+debtadd
    	    			borrowanteG=borrowanteG+1.0
    					borrowG(age,history(age),transhis(age),shockhis(age))=borrowG(age,history(age),transhis(age),shockhis(age)) + 1.0
    					avborintG=avborintG-(1-q(age,history(age),scorehis(age),kprime))*(kprime-bmax-1)/gridfine
    					interestageG(age)=interestageG(age)-(1-q(age,history(age),scorehis(age),kprime))*(kprime-bmax-1)/gridfine
                        profitsG=profitsG-debtadd
                        profitageG(age)=profitageG(age)-debtadd
!                        write (26,998) age,scorehis(age),debtadd,q(age,history(age),scorehis(age),kprime),q(age,history(age),1,kprime)
                    else
    				    debtR=debtR+debtadd
                        debtbyageR(age)=debtbyageR(age)+debtadd
        				borrowanteR=borrowanteR+1.0
    					borrowR(age,history(age),transhis(age),shockhis(age))=borrowR(age,history(age),transhis(age),shockhis(age)) + 1.0
    					avborintR=avborintR-(1-q(age,history(age),scorehis(age),kprime))*(kprime-bmax-1)/gridfine
    					interestageR(age)=interestageR(age)-(1-q(age,history(age),scorehis(age),kprime))*(kprime-bmax-1)/gridfine
                        profitsR=profitsR-debtadd
                        profitageR(age)=profitageR(age)-debtadd
!                        write (27,998) age,scorehis(age),debtadd,q(age,history(age),scorehis(age),kprime),q(age,history(age),scorepoints,kprime)
                    end if
				end if
			else
				debtadd=((kprime-1-bmax)/gridfine)/(1+rs)
				saving=saving+debtadd
				savebyage(age)=savebyage(age)+debtadd
			end if

		end do


		age=lifespan
		history(age)=history(age-1)
		transhis(age)=avtrans
		shockhis(age)=shocks+1
        scorehis(age)=scorehis(age-1)
		wage=phi(history(age))*hbar(age)

		ystig1=(stigma(age-1).eq.1)
		ystig2=(stigma(age-1).eq.2)
		ystig3=(stigma(age-1).eq.3)
		if (ystig1) then
			kpos(age)=bmax+1
            welfarecum=welfarecum + betacum(age-1)*(Utility(autlast+Con1(age,history(age),transhis(age),scorehis(age),shockhis(age)),cigma)-coststig2)
		else if (ystig2) then
			kpos(age)=BRafter(age-1,history(age-1),transhis(age-1),scorehis(age-1),shockhis(age-1))
            welfarecum=welfarecum + betacum(age-1)*Utility(autlast+Con(age,history(age),transhis(age),scorehis(age),kpos(age)),cigma)
		else if (ystig3) then
			kpos(age)=DebtRolled(age-1,history(age-1),transhis(age-1),shockhis(age-1))
            welfarecum=welfarecum + betacum(age-1)*(Utility(autlast+ConB(age,history(age),transhis(age)),cigma)-coststig1)
		else
			kpos(age)=BR(age-1,history(age-1),transhis(age-1),scorehis(age-1),kpos(age-1))
            welfarecum=welfarecum + betacum(age-1)*Utility(autlast+Con(age,history(age),transhis(age),scorehis(age),kpos(age)),cigma)
        end if
        if (kpos(age) < bmax+1) then
			borrowers=borrowers+1.0
			avdoib=avdoib+((bmax+1-kpos(age))/gridfine)/wage
            if (goof) then
        		borrowersG=borrowersG+1.0
    			avdoibG=avdoibG+((bmax+1-kpos(age))/gridfine)/wage
                if (not(ystig3)) then
                    profitsG=profitsG+((bmax+1-kpos(age))/gridfine)/(1+rb)
                    profitageG(age-1)=profitageG(age-1)+((bmax+1-kpos(age))/gridfine)/(1+rb)
                endif
            else
                borrowersR=borrowersR+1.0
            	avdoibR=avdoibR+((bmax+1-kpos(age))/gridfine)/wage
                if (not(ystig3)) then
                    profitsR=profitsR+((bmax+1-kpos(age))/gridfine)/(1+rb)
                    profitageR(age-1)=profitageR(age-1)+((bmax+1-kpos(age))/gridfine)/(1+rb)
                endif
            endif
		end if

		if (wage>0.0) then
			output=output+wage
            if (goof) then
                outputG=outputG+wage
                incoG(age)=incoG(age)+wage
            else
                outputR=outputR+wage
                incoR(age)=incoR(age)+wage
            endif
			varlnincage(age)=varlnincage(age)+(log(wage)-avlnincage(age))**2
		else
			varlnincage(age)=varlnincage(age)+(-1.0d80-avlnincage(age))**2
		end if

		!	Note:	The only bankruptcies here are last periods' defaults		!
		if (ystig3) then
			bankruptcy=bankruptcy+1.0
            if (flag) then
                recidivist=recidivist+1.0
            end if
            if (wage>medianinc) then
                richfiler = richfiler+1.0
            end if
			newbankruptcy(age)=newbankruptcy(age)+1.0
			cons(age)=cons(age)+ConB(age,history(age),transhis(age))+autlast
			avdoid=avdoid+((bmax+1-kpos(age))/gridfine)/wage
			baddebt=baddebt+(bmax+1-kpos(age))/gridfine
			varlnconage(age)=varlnconage(age)+(log(ConB(age,history(age),transhis(age))+autlast)-avlnconage(age))**2
			outdef=outdef+wage
            if (goof) then
    			varlnconG(age)=varlnconG(age)+(log(ConB(age,history(age),transhis(age))+autlast)-avglnconG(age))**2
        		bankruptcyG=bankruptcyG+1.0
				newbankruptcyG(age)=newbankruptcyG(age)+1.0
    			avdoidG=avdoidG+((bmax+1-kpos(age))/gridfine)/wage
				outdefG=outdefG+wage
            else
    			varlnconR(age)=varlnconR(age)+(log(ConB(age,history(age),transhis(age))+autlast)-avglnconR(age))**2
        		bankruptcyR=bankruptcyR+1.0
				newbankruptcyR(age)=newbankruptcyR(age)+1.0
    			avdoidR=avdoidR+((bmax+1-kpos(age))/gridfine)/wage
				outdefR=outdefR+wage
            endif
			increason=1		!	No income shock	!
			transreason=1		!	No income shock	!
			BRbyreason(shockhis(age),increason,transreason)=BRbyreason(shockhis(age),increason,transreason)+1.0
            if (goof) then
                BRbyreasonG(shockhis(age),increason,transreason)=BRbyreasonG(shockhis(age),increason,transreason)+1.0
            else
                BRbyreasonR(shockhis(age),increason,transreason)=BRbyreasonR(shockhis(age),increason,transreason)+1.0
            end if
		else if (ystig1) then
			cons(age)=cons(age)+Con1(age,history(age),transhis(age),scorehis(age),shockhis(age))+autlast
			varlnconage(age)=varlnconage(age)+(log(Con1(age,history(age),transhis(age),scorehis(age),shockhis(age))+autlast)-avlnconage(age))**2
            if (goof) then
    			varlnconG(age)=varlnconG(age)+(log(Con1(age,history(age),transhis(age),scorehis(age),shockhis(age))+autlast)-avglnconG(age))**2
            else
    			varlnconR(age)=varlnconR(age)+(log(Con1(age,history(age),transhis(age),scorehis(age),shockhis(age))+autlast)-avglnconR(age))**2
            end if
		else
			cons(age)=cons(age)+Con(age,history(age),transhis(age),scorehis(age),kpos(age))+autlast
			varlnconage(age)=varlnconage(age)+(log(Con(age,history(age),transhis(age),scorehis(age),kpos(age))+autlast)-avlnconage(age))**2
            if (goof) then
    			varlnconG(age)=varlnconG(age)+(log(Con(age,history(age),transhis(age),scorehis(age),kpos(age))+autlast)-avglnconG(age))**2
            else
    			varlnconR(age)=varlnconR(age)+(log(Con(age,history(age),transhis(age),scorehis(age),kpos(age))+autlast)-avglnconR(age))**2
            end if
		endif
        realizedEV = realizedEV + welfarecum
        if (goof) then
            realizedG=realizedG+welfarecum
        else
            realizedR=realizedR+welfarecum
        end if
        if (welfarecum>bestluck) then
            bestluck = welfarecum
        end if
        if (welfarecum<worstluck) then
            worstluck = welfarecum
        end if

	end do

    realizedwelfare(j) = realizedEV/Pop
	bankrtotal(j)=(bankruptcy/Pop)/lifespan
    netprofitG(j)=(profitsG/Pop)/lifespan
    netprofitR(j)=(profitsR/Pop)/lifespan
    if (countR>0.0) then
        welfareR(j)=realizedR/countR
        bankrR(j)=(bankruptcyR/countR)/lifespan
    end if
    if (countG>0.0) then
        welfareG(j)=realizedG/countG
        bankrG(j)=(bankruptcyG/countG)/lifespan
        overborrowtotal(j)=(overborrowing/countG)/lifespan
        overborrqtotal(j)=(overborrq/countG)/lifespan
        overhangtotal(j)=(overhang/countG)/lifespan
        gaintotal(j)=(gain/countG)/lifespan
        totalover(j)=(countover/countG)/(lifespan-1)
        totalunder(j)=(countunder/countG)/(lifespan-1)
        totalearly(j)=(tooearly/countG)/(lifespan-1)
        totallate(j)=(toolate/countG)/(lifespan-1)
    end if
	defaulttotal(j)=(defaults/Pop)/lifespan
	if (debt>0.0) then
		avinttotal(j)=avborint/debt
		debttotal(j)=(debt/Pop)/lifespan
		borrowerantetotal(j)=(borrowante/Pop)/lifespan
	else
		avinttotal(j)=0.0
		debttotal(j)=0.0
		borrowerantetotal(j)=0.0
	end if
	if (debtG>0.0) then
		avintG(j)=avborintG/debtG
		debttotalG(j)=(debtG/countG)/lifespan
		borroweranteG(j)=(borrowanteG/countG)/lifespan
	else
		avintG(j)=0.0
		debttotalG(j)=0.0
		borroweranteG(j)=0.0
	end if
	if (debtR>0.0) then
		avintR(j)=avborintR/debtR
		debttotalR(j)=(debtR/countR)/lifespan
		borroweranteR(j)=(borrowanteR/countR)/lifespan
	else
		avintR(j)=0.0
		debttotalR(j)=0.0
		borroweranteR(j)=0.0
	end if
	savetotal(j)=(saving/Pop)/lifespan
	if (borrowers>0.0) then
		borrowerposttotal(j)=(borrowers/Pop)/lifespan
		debtincbortotal(j)=avdoib/borrowers
	else
		borrowerposttotal(j)=0.0
		debtincbortotal(j)=0.0
	end if
	if (borrowersG>0.0) then
		borrowerpostG(j)=(borrowersG/countG)/lifespan
		debtincborG(j)=avdoibG/borrowersG
	else
		borrowerpostG(j)=0.0
		debtincborG(j)=0.0
	end if
	if (borrowersR>0.0) then
		borrowerpostR(j)=(borrowersR/countR)/lifespan
		debtincborR(j)=avdoibR/borrowersR
	else
		borrowerpostR(j)=0.0
		debtincborR(j)=0.0
	end if
	if (bankruptcy>0.0) then
		debtincdeftotal(j)=avdoid/bankruptcy
		writeoff(j)=baddebt/bankruptcy
		savers(j)=wannasave/bankruptcy
        repeatfrac(j)=recidivist/bankruptcy
        richfrac(j)=richfiler/bankruptcy
	else
		debtincdeftotal(j)=0.0
		writeoff(j)=0.0
		savers(j)=0.0
        repeatfrac(j)=0.0
        richfrac(j)=0.0
	end if
	if (bankruptcyG>0.0) then
		debtincdefG(j)=avdoidG/bankruptcyG
	else
		debtincdefG(j)=0.0
	end if
	if (bankruptcyR>0.0) then
		debtincdefR(j)=avdoidR/bankruptcyR
	else
		debtincdefR(j)=0.0
	end if

end do
!close (26)
!close (27)

varlnconage = (varlnconage/Pop)/(Monte-Monte1)
varlnconG = (varlnconG/Pop)/(Monte-Monte1)
varlnconR = (varlnconR/Pop)/(Monte-Monte1)
varlnincage = (varlnincage/Pop)/(Monte-Monte1)
profitageG = (profitageG/Pop)/Monte
profitageR = (profitageR/Pop)/Monte
varlncon=0.0
varlninc=0.0
do age=1,lifespan
	varlncon = varlncon + varlnconage(age)
	varlninc = varlninc + varlnincage(age)
end do
varlncon = varlncon/lifespan
varlninc = varlninc/lifespan

debt=0.0
do i=(Monte1+1),Monte
    debt = debt + debttotal(i)
end do
varint = (varint/debt)/(Pop*lifespan*(Monte-Monte1))

994 format (I4,14F9.5)
open (18,file='MonteOut.txt')
write (18,*) 'Input ID#',id
write (18,*) Monte,'experiments with',Pop,'people in each generation'
write (18,*) 'Run Bankruptcies AvInterest Debt Savings ExAnteBorrowers ExPostBorrowers (Debt/Income)Borrowers (Debt/Income)Bankrupts BadDebt RealizedWelfare OverBorrowing(flow.face) FileTooLate GoofyProfit RealistProfit'
do i=1,Monte
write (18,994) i, bankrtotal(i), avinttotal(i), debttotal(i), savetotal(i), borrowerantetotal(i), borrowerposttotal(i), debtincbortotal(i), debtincdeftotal(i), writeoff(i), realizedwelfare(i), overborrowtotal(i), totallate(i), netprofitG(i), netprofitR(i)
end do
close (18)

995 format (I4,8F9.5)
open (17,file='MonteGoofy.txt')
write (17,*) 'Input ID#',id
write (17,*) Monte,'experiments with',Pop,'people in each generation'
write (17,*) 'Results for the goofy portion of the population, whose share is', lambda
write (17,*) 'Run BankruptRate AvInterest Debt ExAnteBorrowers ExPostBorrowers (Debt/Income)Borrowers (Debt/Income)Bankrupts RealizedWelfare'
do i=1,Monte
	write (17,995) i, bankrG(i), avintG(i), debttotalG(i), borroweranteG(i), borrowerpostG(i), debtincborG(i), debtincdefG(i), welfareG(i)
end do
close (17)

open (16,file='MonteRealist.txt')
write (16,*) 'Input ID#',id
write (16,*) Monte,'experiments with',Pop,'people in each generation'
write (16,*) 'Results for the realist portion of the population'
write (16,*) 'Run BankruptRate AvInterest Debt ExAnteBorrowers ExPostBorrowers (Debt/Income)Borrowers (Debt/Income)Bankrupts RealizedWelfare'
do i=1,Monte
	write (16,995) i, bankrR(i), avintR(i), debttotalR(i), borroweranteR(i), borrowerpostR(i), debtincborR(i), debtincdefR(i), welfareR(i)
end do
close (16)


debt=0.0
debtG=0.0
debtR=0.0
saving=0.0
bankruptcy=0.0
recidivist=0.0
richfiler=0.0
bankruptcyG=0.0
bankruptcyR=0.0
wannasave=0.0
defaults=0.0
avborint=0.0
avborintG=0.0
avborintR=0.0
borrowante=0.0
borrowanteG=0.0
borrowanteR=0.0
borrowers=0.0
borrowersG=0.0
borrowersR=0.0
avdoid=0.0
avdoidG=0.0
avdoidR=0.0
avdoib=0.0
avdoibG=0.0
avdoibR=0.0
baddebt=0.0
realizedEV=0.0
realizedG=0.0
realizedR=0.0
overborrowing=0.0
overborrq=0.0
overhang=0.0
gain=0.0
countover=0.0
countunder=0.0
tooearly=0.0
toolate=0.0
profitsG=0.0
profitsR=0.0
do i=1,Monte
	debt=debt+debttotal(i)
	debtG=debtG+debttotalG(i)
	debtR=debtR+debttotalR(i)
	saving=saving+savetotal(i)
	avborint=avborint+avinttotal(i)*debttotal(i)
	avborintG=avborintG+avintG(i)*debttotalG(i)
	avborintR=avborintR+avintR(i)*debttotalR(i)
	bankruptcy=bankruptcy+bankrtotal(i)
    recidivist=recidivist+repeatfrac(i)
    richfiler=richfiler+richfrac(i)
	bankruptcyG=bankruptcyG+bankrG(i)
	bankruptcyR=bankruptcyR+bankrR(i)
	defaults=defaults+defaulttotal(i)
	borrowante=borrowante+borrowerantetotal(i)
	borrowanteG=borrowanteG+borroweranteG(i)
	borrowanteR=borrowanteR+borroweranteR(i)
	borrowers=borrowers+borrowerposttotal(i)
	borrowersG=borrowersG+borrowerpostG(i)
	borrowersR=borrowersR+borrowerpostR(i)
	avdoid=avdoid+debtincdeftotal(i)*bankrtotal(i)
	avdoidG=avdoidG+debtincdefG(i)*bankrG(i)
	avdoidR=avdoidR+debtincdefR(i)*bankrR(i)
	baddebt=baddebt+writeoff(i)*bankrtotal(i)
	wannasave=wannasave+savers(i)*bankrtotal(i)
	avdoib=avdoib+debtincbortotal(i)*borrowerposttotal(i)
	avdoibG=avdoibG+debtincborG(i)*borrowerpostG(i)
	avdoibR=avdoibR+debtincborR(i)*borrowerpostR(i)
    realizedEV=realizedEV+realizedwelfare(i)
    realizedG=realizedG+welfareG(i)
    realizedR=realizedR+welfareR(i)
    overborrowing=overborrowing+overborrowtotal(i)
    overborrq=overborrq+overborrqtotal(i)
    overhang=overhang+overhangtotal(i)
    gain=gain+gaintotal(i)
    countover=countover+totalover(i)
    countunder=countunder+totalunder(i)
    tooearly=tooearly+totalearly(i)
    toolate=toolate+totallate(i)
    profitsG=profitsG+netprofitG(i)
    profitsR=profitsR+netprofitR(i)
end do
avborint=avborint/debt		!	Debt-weighted across simulations!
if (debtG>0.0) then
    avborintG=avborintG/debtG
else
    avborintG=0.0
endif
if (debtR>0.0) then
    avborintR=avborintR/debtR
else
    avborintR=0.0
endif

debt=debt/Monte
debtG=debtG/Monte
debtR=debtR/Monte
saving=saving/Monte
bankruptcy=bankruptcy/Monte
recidivist=recidivist/Monte
richfiler=richfiler/Monte
bankruptcyG=bankruptcyG/Monte
bankruptcyR=bankruptcyR/Monte
defaults=defaults/Monte
borrowante=borrowante/Monte
borrowanteG=borrowanteG/Monte
borrowanteR=borrowanteR/Monte
borrowers=borrowers/Monte
borrowersG=borrowersG/Monte
borrowersR=borrowersR/Monte
realizedEV=realizedEV/Monte
realizedG=realizedG/Monte
realizedR=realizedR/Monte
overborrowing=overborrowing/Monte
overborrq=overborrq/Monte
overhang=overhang/Monte
gain=gain/Monte
countover=countover/Monte
countunder=countunder/Monte
tooearly=tooearly/Monte
toolate=toolate/Monte
profitsG=profitsG/Monte
profitsR=profitsR/Monte

cons=(cons/Pop)/Monte
debtbyage=(debtbyage/Pop)/Monte
do age=1,(lifespan-1)
    if (debtbyageG(age)>0.0) then
        interestageG(age)=interestageG(age)/debtbyageG(age)
    else
        interestageG(age)=0.0
    end if
    if (debtbyageR(age)>0.0) then
        interestageR(age)=interestageR(age)/debtbyageR(age)
    else
        interestageR(age)=0.0
    end if
end do
!
!   Rather than dividing by Pop*Monte*lambda below, we should really just divide by PopG
!
!                       FIX WHEN YOU GET A CHANCE
!
if (lambda>0.0) then
    debtbyageG=(debtbyageG/(Pop*lambda))/Monte
end if
debtbyageR=(debtbyageR/(Pop*(1.0-lambda)))/Monte
savebyage=(savebyage/Pop)/Monte
newbankruptcy=(newbankruptcy/Pop)/Monte
newbankruptcyG=(newbankruptcyG/(Pop*lambda))/Monte
newbankruptcyR=(newbankruptcyR/(Pop*(1-lambda)))/Monte
toolatebyage=(toolatebyage/Pop)/Monte
avgscoreG=avgscoreG/(PopG*scorebins)
avgscoreR=avgscoreR/(PopR*scorebins)
BRbyreason=(BRbyreason/Pop)/(Monte*lifespan)
BRbyreasonG=(BRbyreasonG/(Pop*lambda))/(Monte*lifespan)
BRbyreasonR=(BRbyreasonR/(Pop*(1-lambda)))/(Monte*lifespan)

output=(output/Pop)/(Monte*lifespan)
outputG=(outputG/(Pop*lambda))/(Monte*lifespan)
outputR=(outputR/(Pop*(1-lambda)))/(Monte*lifespan)
incoG=(incoG/(Pop*lambda))/Monte
incoR=(incoR/(Pop*(1-lambda)))/Monte
if (bankruptcy>0) then
	outdef = (outdef/Pop)/(bankruptcy*Monte*lifespan)
	avdoid = avdoid/(bankruptcy*Monte)
	baddebt = baddebt/(bankruptcy*Monte)
	wannasave = wannasave/(bankruptcy*Monte)
end if
if (bankruptcyG>0.0) then
	outdefG = (outdefG/(Pop*lambda))/(bankruptcyG*Monte*lifespan)
	avdoidG = avdoidG/(bankruptcyG*Monte)
end if
if (bankruptcyR>0.0) then
	outdefR = (outdefR/(Pop*(1-lambda)))/(bankruptcyR*Monte*lifespan)
	avdoidR = avdoidR/(bankruptcyR*Monte)
end if
if (borrowers>0) then
	avdoib = avdoib/(borrowers*Monte)
end if
if (borrowersG>0.0) then
	avdoibG = avdoibG/(borrowersG*Monte)
end if
if (borrowersR>0.0) then
	avdoibR = avdoibR/(borrowersR*Monte)
end if

do age=1,(lifespan-1)
	do s=1,states
		do st=1,trans
			do sh=1,shocks+1
				if (borrow(age,s,st,sh)>0.0) then
					borrowbyage(age)=borrowbyage(age)+borrow(age,s,st,sh)
					borrowbyageG(age)=borrowbyageG(age)+borrowG(age,s,st,sh)
					borrowbyageR(age)=borrowbyageR(age)+borrowR(age,s,st,sh)
!					asum=prob0(s)*probtTrue(st)*probsTrue(sh)*Pop*Monte
!					fracborrow(age,s,st,sh)=borrow(age,s,st,sh)/asum
				end if
			end do
		end do
	end do
	borrowbyage(age)=(borrowbyage(age)/Pop)/Monte
	borrowbyageG(age)=(borrowbyageG(age)/(Pop*lambda))/Monte
	borrowbyageR(age)=(borrowbyageR(age)/(Pop*(1-lambda)))/Monte
end do

do i=1,groups
	debincB(i)=debtB(i)/outputB(i)
	outputB(i)=outputB(i)/PopB(i)
	debtB(i)=debtB(i)/PopB(i)
	borrowB(i)=borrowB(i)/PopB(i)
	avdebincB(i)=avdebincB(i)/borpostB(i)
	debtpostB(i)=debtpostB(i)/PopB(i)
	borpostB(i)=borpostB(i)/PopB(i)
	filingsB(i)=filingsB(i)/PopB(i)
end do


open (20,file='outputlong.txt')
write (20,*) 'Input ID#',id
write (20,*) 'Exemption, Garnishment, Burnt:', cbar,garn,burn
write (20,*) 'Autonomous, Saving rate, Risk-free borrowing rate, Roll-over rate:', aut,rs,rb,rroll
write (20,*) 'Expense shocks intercepts:',expenconst
write (20,*) 'Expense shocks slopes:',expenprop
write (20,*) 'Results:'
write (20,*) 'Perceived EV=',EV
write (20,*) 'Paternalistic Goofy EV=',EVpat
write (20,*) 'Educated Goofy EV=',EVedu
write (20,*) 'Realized Average Welfare=',realizedEV
write (20,*) 'Realized Goofy Welfare=',realizedG
write (20,*) 'Realized Realist Welfare=',realizedR
write (20,*) 'Luckiest realized utility=',bestluck
write (20,*) 'Unluckiest realized utility=',worstluck
write (20,*) 'Fraction declaring bankruptcy',bankruptcy
write (20,*) 'Fraction of filings from repeat filers',recidivist
write (20,*) 'Fraction of filers with income above the median',richfiler
write (20,*) 'Bankruptcy rate among the goofy',bankruptcyG
write (20,*) 'Bankruptcy rate among the realists',bankruptcyR
write (20,*) 'Fraction in default (after bankruptcy)',defaults
write (20,*) 'Fraction of goofy who would have filed if aware',toolate
write (20,*) 'Fraction of goofy who filed, but would NOT have if aware',tooearly
write (20,*) 'Total earnings =',output
write (20,*) 'Goofy earnings =',outputG
write (20,*) 'Debt=',debt
write (20,*) 'Goofy Debt=',debtG
write (20,*) 'Realist Debt=',debtR
write (20,*) 'Overborrowing (avg flow in face value)=',overborrowing
write (20,*) 'Overborrowing (flow with prices)=',overborrq
write (20,*) 'Overborrowing from not filing (with prices)=',overhang
write (20,*) 'Fraction of goofy overborrowing =',countover
write (20,*) 'Fraction of goofy underborrowing =',countunder
write (20,*) 'Welfare gain from education (in utils)=',gain
write (20,*) 'Net profits from Goofy=',profitsG
write (20,*) 'Net profits from Realists=',profitsR
write (20,*) 'Debt/Earnings=',debt/output
write (20,*) 'Goofy Debt/Earnings=',debtG/outputG
write (20,*) 'Realist Debt/Earnings=',debtR/outputR
write (20,*) 'Savings=',saving
write (20,*) 'Savings/Earnings=',saving/output
write (20,*) 'Fraction of bankrupts who want to save',wannasave
write (20,*) 'Fraction of ex-post borrowers = ',borrowers
write (20,*) 'Fraction of ex-post borrowers among goofy = ',borrowersG
write (20,*) 'Fraction of ex-post borrowers among realists = ',borrowersR
write (20,*) 'Fraction of ex-ante borrowers = ',borrowante
write (20,*) 'Fraction of ex-ante borrowers among goofy = ',borrowanteG
write (20,*) 'Fraction of ex-ante borrowers among realists = ',borrowanteR
if (bankruptcy>0) then
	write (20,*) 'Mean Gross Earnings of Defaulters=',outdef
	write (20,*) 'BankruptY / AverageY =',outdef/output
	write (20,*) 'Average (ex-post) debt/income of bankrupts:',avdoid
	write (20,*) 'Ex-post debt of bankrupts:',baddebt
	write (20,*) 'Av(ex-post)debt/Av.income of bankrupts:',baddebt/outdef
end if
if (bankruptcyG>0.0) then
	write (20,*) 'Mean Gross Earnings of Goofy Defaulters=',outdefG
	write (20,*) 'Average (ex-post) debt/income of Goofy bankrupts:',avdoidG
end if
if (bankruptcyR>0.0) then
	write (20,*) 'Mean Gross Earnings of Realist Defaulters=',outdefR
	write (20,*) 'Average (ex-post) debt/income of Realist bankrupts:',avdoidR
end if
write (20,*) 'Average (ex-post) debt/income of borrowers:',avdoib
write (20,*) 'Average (ex-post) debt/income of goofy borrowers:',avdoibG
write (20,*) 'Average (ex-post) debt/income of realist borrowers:',avdoibR
write (20,*) 'Average borrowing interest rate:',avborint
write (20,*) 'Average borrowing interest rate of goofy:',avborintG
write (20,*) 'Average borrowing interest rate of realists:',avborintR
write (20,*) 'Variance of borrowing interest rates:',varint
write (20,*) 'Coefficient of variation of borrowing interest rates:',sqrt(varint)/avborint
write (20,*) 'Variance of log Consumption (incl. autonomous):',varlncon
write (20,*) 'Variance of log Income:',varlninc

write (20,*) 'Profits from goofy by age:',profitageG
write (20,*) 'Profits from realists by age:',profitageR
write (20,*) 'New Bankruptcies by age:',newbankruptcy
write (20,*) 'Debt (ex-ante) by age:',debtbyage
write (20,*) 'Savings (before interest) by age:',savebyage
write (20,*) 'Fraction of borrowers (ex-ante) by age:',borrowbyage
write (20,*) 'Average Consumption (incl. autonomous) by age:',cons
write (20,*) 'Variance of log Consumption (incl. autonomous) by age:',varlnconage
write (20,*) 'Variance of log Consumption of goofy:',varlnconG
write (20,*) 'Variance of log Consumption of realists:',varlnconR
write (20,*) 'Variance of log Income by age:',varlnincage

write (20,*) 'Bankruptcies by reason:'
write (20,*) 'No bad Markov shock, No bad transitory shock:',BRbyreason(1:(shocks+1),1,1)
write (20,*) 'Bad Markov shock, No bad transitory shock:',BRbyreason(1:(shocks+1),2,1)
write (20,*) 'No bad Markov shock, Bad transitory shock:',BRbyreason(1:(shocks+1),1,2)
write (20,*) 'Bad Markov shock, Bad transitory shock:',BRbyreason(1:(shocks+1),2,2)
write (20,*) 'Bankruptcies by reason among Goofy:'
write (20,*) 'No bad Markov shock, No bad transitory shock:',BRbyreasonG(1:(shocks+1),1,1)
write (20,*) 'Bad Markov shock, No bad transitory shock:',BRbyreasonG(1:(shocks+1),2,1)
write (20,*) 'No bad Markov shock, Bad transitory shock:',BRbyreasonG(1:(shocks+1),1,2)
write (20,*) 'Bad Markov shock, Bad transitory shock:',BRbyreasonG(1:(shocks+1),2,2)
write (20,*) 'Bankruptcies by reason among Realists:'
write (20,*) 'No bad Markov shock, No bad transitory shock:',BRbyreasonR(1:(shocks+1),1,1)
write (20,*) 'Bad Markov shock, No bad transitory shock:',BRbyreasonR(1:(shocks+1),2,1)
write (20,*) 'No bad Markov shock, Bad transitory shock:',BRbyreasonR(1:(shocks+1),1,2)
write (20,*) 'Bad Markov shock, Bad transitory shock:',BRbyreasonR(1:(shocks+1),2,2)
close (20)


open (22,file='NewOrderPlus.txt')
write (22,*) id
write (22,*) realizedR
write (22,*) realizedG
write (22,*) realizedEV
write (22,*) bankruptcyR
write (22,*) bankruptcyG
write (22,*) bankruptcy
write (22,*) toolate
write (22,*) tooearly
write (22,*) avborintR
write (22,*) avborintG
write (22,*) avborint
write (22,*) debtR/outputR
write (22,*) debtG/outputG
write (22,*) debt/output
write (22,*) overborrowing
write (22,*) avdoibR
write (22,*) avdoibG
write (22,*) avdoib
if (bankruptcyR>0.0) then
	write (22,*) avdoidR
else
    write (22,*) 'No realist bankruptcy'
end if
if (bankruptcyG>0.0) then
	write (22,*) avdoidG
else
    write (22,*) 'No goofy bankruptcy'
end if
if (bankruptcy>0) then
	write (22,*) avdoid
else
    write (22,*) 'No bankruptcy'
end if
write (22,*) outputR
write (22,*) outputG
write (22,*) output
write (22,*) debtR
write (22,*) debtG
write (22,*) debt
write (22,*) borrowersR
write (22,*) borrowersG
write (22,*) borrowers
write (22,*) borrowanteR
write (22,*) borrowanteG
write (22,*) borrowante
if (bankruptcyR>0.0) then
	write (22,*) outdefR
else
    write (22,*) 'No realist bankruptcy'
end if
if (bankruptcyG>0.0) then
	write (22,*) outdefG
else
    write (22,*) 'No goofy bankruptcy'
end if
if (bankruptcy>0) then
	write (22,*) outdef
else
    write (22,*) 'No bankruptcy'
end if
write (22,*) saving
write (22,*) defaults
write (22,*) outdef/output
write (22,*) baddebt
write (22,*) baddebt/outdef
write (22,*) varlncon
write (22,*) varlninc
write (22,*) profitsR
write (22,*) profitsG
write (22,*) recidivist
write (22,*) sqrt(varint)/avborint
write (22,*) gain
write (22,*) overborrq
write (22,*) overhang
write (22,*) countover
write (22,*) countunder
write (22,*) EV
write (22,*) EVpat
write (22,*) EVedu
close (22)


997 format (I3,15F9.5,2F10.7)
open (19,file='Life-Cycle.txt')
write (19,*) 'Input ID#',id
write (19,*) ' Age (Ex-ante)DebtGoofy (Ex-ante)DebtRealist IncomeGoofy IncomeRealist (Ex-ante)BorrowersGoofy (Ex-ante)BorrowersRealist BankruptcyRateGoofy BankruptcyRateRealist AvgInterestGoofy AvgInterestRealist AvgScoreGoofy AvgScoreRealist WouldHaveFiled VarLnConGoofy VarLnConRealist NetProfitGoofy NetProfitRealist'
do age=1,(lifespan-1)
    write (19,997) age, debtbyageG(age), debtbyageR(age), incoG(age), incoR(age), borrowbyageG(age), borrowbyageR(age), newbankruptcyG(age), newbankruptcyR(age), interestageG(age), interestageR(age), avgscoreG(age), avgscoreR(age), toolatebyage(age), varlnconG(age), varlnconR(age), profitageG(age), profitageR(age)
end do
close (19)


993 format (I4,9F9.5)
open (23,file='BinsOut.txt')
write (23,*) 'Input ID#',id
write (23,*) ' Bin Income (Ex-ante)Debt (Ex-ante)Borrowers (Ex-ante)Debt/Income (Ex-post)Debt (Ex-post)Borrowers (Ex-post)Debt/Income Av(Ex-post)Debt/Income_of_borrowers FractionFiling'
do i=1,groups
	write (23,993) i, outputB(i), debtB(i), borrowB(i), debincB(i), debtpostB(i), borpostB(i), (debtpostB(i)/outputB(i)), avdebincB(i), filingsB(i)
end do
close (23)

!wss=PopG+PopR
!print *,'Fraction if goofy',PopG/wss
print *,'Fraction declaring bankruptcy',bankruptcy
print *,'Fraction of filings from repeat filers',recidivist
print *,'Fraction of filers with income above the median',richfiler
print *,'Fraction in default (after bankruptcy)',defaults
print *,'New Bankruptcies:',newbankruptcy
print *,'Debt=',debt
print *,'Earnings=',output
print *,'Debt/Earnings=',debt/output
print *,'Savings=',saving
print *,'Savings/Earnings=',saving/output
print *,'Fraction of ex-post borrowers = ',borrowers
print *,'Fraction of ex-ante borrowers = ',borrowante
if (bankruptcy>0.0) then
	print *,'Mean Gross Earnings of Defaulters=',outdef
	print *,'BankruptY / AverageY =',outdef/output
	print *,'Average (ex-post) debt/income of bankrupts =',avdoid
	print *,'Ex-post debt of bankrupts:',baddebt
	print *,'Av(ex-post)debt/Av.income of bankrupts:',baddebt/outdef
end if
print *,'Average (ex-post) debt/income of borrowers =',avdoib
print *,'Realists EV=',EV
!print *,'Aware (paternalistic) Goofy EV=',EVedu
print *,'Fraction of goofy who would have filed if aware',toolate
print *,'Fraction of goofy who filed, but would NOT have if aware',tooearly

!pause
777 end program



Real(8) Function Utility(c,cig)
real(8) c, cig
Utility=c**(1-cig)/(1-cig)
end function
