' This program estimates stochastic VAR regressions as described in Ferrucci and Penalver (2003).
' Code written by Marcos Rietti Souto, while he worked at the IMF Research Department.
' It is set to work with quarterly data. With few modifications, it can work with different data frequencies.

' This version: August 9, 2005.

' Inputting the data.

cd "C:\Documents and Settings\xdebrun\My Local Documents\RESEARCH\CDO(2005)\Simulations"			' This line needs to be modified, depending on your system.

db temp

!Firstobs= 1								' Sets the initial obs of the TS, for setting up the workfile.

!Lastobs=100							' Sets the last obs of the TS, for setting up the workfile.

workfile xvar u !Firstobs !Lastobs

Call Setup

Call GetX

Call GetLastQ

Call HP_TS

Call GetVAR

Call CalcInitialCond

Call StochVAR

Call Annualize

Call ForGapPSDebt

Call GetTrend

Call GetProb

Call CalcStat


'****************************************************************************************************************************************
' This subroutine sets up parameters for the entire code. It needs to be changed according to the data being used.
Subroutine Setup

%Country_Name = "Countryname"					' Set here the name of the country you are working with.
 
%x_variables = "obs year quarter GDP int_rate REER commodity USrate"	
														' Variables names to be included in the VAR.
														' First 3 columns are always the same (obs., year, and quarter).
														' The names have to be consistent with the !nTS (number of Time Series)
														' parameter below. ALWAYS name real gdp level time series as GDP,
														' otherwise HP filtering below won't work.

%x_var = "USrate int_rate growth l_REER"	
														' Here put the variables in the sequence you want them to be estimated
														' in the VAR, from the most exogeneous to the most endogeneous.
														' This sequence matters for the Cholesky factorization. Now, put growth
														' in place of GDP (VAR will be estimated using GDP growth and not
														' GDP levels.
														' %x_var also needs to be consistent with the vector annualize below
														' (only data in %x_var will be estimated and later annualized - see
														' annualize vector below for explanation on how to set it).
														' %x_var can be a sub-set of %x_variables and NOT include ALL
														' variables that are in %x_variables.
														' IF WE CHANGE THIS SEQUENCE, WE HAVE TO CHANGE DEBT
														' EXPRESSION ACCORDINGLY (INDEX REFERENCES).

' Input data must be a in a txt file in the following format. It needs to have a headline, which will ne skipped when reading
' the file below. If original data does not have a headline, then needs to chnage the reading command below.
' Obs#		Year	Quarter	Var1			Var2		Var3	...
' 1			1994	1			32000			0.070
' 2			1994	2			33361			0.068
' 3			1994	3			33011			0.071
' 4			1994	4			31933			0.085
' 5			1995	1			32255			0.076


!nTS = 4									' Number of variables in the x vector. Note that Obs#, Year, and Quarter DO NOT
											' count for !nTS

vector(!nTS) annualize					' This vector sets the choice for annualization for each of the variables above:
annualize.fill  1, 1, 1, 0				' 1 is for compounding and 0 is for averaging. It MUST have the same numbers
											' of 0's and 1's, as the number of variables in the input dataset (!nTS).

!debt_choice = 1							' If debt_choice = 0 then use simple debt expression, otherwise, use the more
											' complete version. If using more complete version, fraction of foreign debt over
											' forecast periods needs to be specified below.

if !debt_choice = 0 then
	%x_var = "int_rate growth l_REER"
	!nTS = 3
	vector(!nTS) annualize
	annualize.fill  1, 1, 0
endif


!NHP=1600								' Parameter for HP estimation. !NHP = 100 for annual, 1600 for 
											' quarterly, and 14400 for monthly.

!Maxlags = 4								' Maximum number of lags to be used in the optmizing Akaike criterion.

!Ch_lag = 0								' If ch_lag = 1 then use automatic procedure, otherwise, use the !setlag below

!setlag = 1								' Number of lags in the VAR if do not want to use minimizing Akaike criterion

!Nfor = 20									' Number of quarterly forecast periods to be used in the stochastic simulation.

!Nfor_a = @ceiling(!Nfor/4)				' Approximate number of annual forecast periods (for PS and debt forecasts).
											' It is said approximate because if !Nfor is not a multiple of 4, then
											' @ceiling(!Nfor/4) will be different than !Nfor/4. @ceiling(!Nfor/4) will get the
											' FIRST integer !Nfor_a, such that !Nfor_a >= @ceiling(!Nfor/4).

!nSim = 1000								' Number of simulation runs.

!IniCond = 1								' Will determine the initial condition for stochastic VAR: 1 for steady-state and 2 
											' for last historical observation. If there are missing observations for any of the
											' historical lags used for forecasting, then the program automatically shuts
											' the option 2 off, and uses the steady-state condition.

!outp = 1									' If 1, output debt forecasts to txt files. For any other value, it won't ouput.

!GapChoice = 0							' If 1, code will use the GDP quarterly data to estimate the last gap. For any
											' other value, it will use the !lastgap provided by user below.

!lastgap = 0						' This inserts the last calculated output gap, as a fraction of GDP level, for 
											' forecasting output gaps.

' This is being calculated below, from the steady state initial condition, but I will keep the option of setting this as a
' parameter (I just have it commented out)
'!potgrowth = 0.035						' This inserts the potential growth for the GDP trend.

!Curdebt = 0.545							' This inserts the current stock of debt, for forecasting debt.

!debtfirstyr = 1995						' First year of historical data on debt.

!PS_Choice = 0							' Allows for inserting different arbitrary shocks in PS equation, in each
											' forecasting period. If !PS_Choice = 0 then no arbitrary shock is inserted
											' and if !PS_Choice = 0 then includes the arbitrary shocks as below.

vector(!Nfor_a) PS_Shocks				' This is the vector with arbitrary shocks to PS equation.

PS_Shocks.fill  -0.01, -0.01, -0.01, -0.01, -0.01			
											' These numbers need to be consistent with !Nfor_a. If !Nfor_a = 5 then we
											' need 5 numbers after annualize. PS_Shocks, each one of them
											' representing the arbitrary shock in each forecast period, separated by comma.

!c_ps = 0.01804									' Constant from PS regression.

!b1_ps =0.32790								' Gap coefficient from PS regression.

!b2_ps =0.04614								' Debt coefficient from PS regression.

!var_ps =0.00029					' Variance of error term from PS regression.

vector(!Nfor_a) debt_s					' Vector with foreign debt fraction for each forecast period. Numbers are provided below.

debt_s.fill  0.2294, 0.2294, 0.2294, 0.2294, 0.2294			


!perc = 0.9								' Estimates the debt threshold below which falls 90% of the simulated debt
											' in each forecasted period.

!upperdebt = 0.645							' For estimating the probability that debt may go higher than this threshold.

!lowerdebt = 0.545							' For estimating the probability that debt may fall below this threshold.

%XSourceData1 = "data\XXX_1.txt"		' This provides the path, in adition to cd above, where the input data is. It will be
											' used in the read statement below.

%XSourceData2 = "data\XXXdebt.txt"

%OutputPath = "data"					' This provides the path to which the output will be saved

Endsub		' End of subroutine Setup.
'****************************************************************************************************************************************



'****************************************************************************************************************************************
' This subroutine inputs the data.
Subroutine GetX

matrix(!Lastobs,!nTS+3) x_data
x_data.read(skiprow=1,t=txt, d=t) {%XSourceData1}		' Reads the X data in matrix format from a txt file. Will be used to 
																' find the last quarter for which there is data for all variables.

Read(t=txt, d=t) {%XSourceData1} {%x_variables}	    	' Reads the data from a txt file to run the VAR. It creates series 
																' for each TS, which is the appropriate format for running the VAR.

Endsub		' End of subroutine GetX.
'****************************************************************************************************************************************



'****************************************************************************************************************************************
' This subroutine gets the last period for which there is data for all variables.
Subroutine GetLastQ

' This is to find the last row in the annual matrix, for which there is data for all TS.
for !n = !Lastobs to 1 step -1
	!count1 = 0
	for !n2 = 1 to !nTS
		if x_data(!n,!n2+3) <> NA then 
			!count1 = !count1 + 1
		endif
	next	' !n2
	if !count1 = !nTS then
		!lasto = x_data(!n,1)
		!lasty = x_data(!n,2)
		!lastq = x_data(!n,3)
		exitloop
	endif
next	' !n

' The !lastq parameter can take the values of 1, 2, 3, or 4, representing 1st, 2nd, 3rd, or 4th quarters. The !lasty
' parameter represents the last year of data for all TS. The "lasto gets the location of the last observation, which
' will be used down below for compounding growth rates.

Endsub		' End of subroutine GetLastQ.
'****************************************************************************************************************************************



'****************************************************************************************************************************************
' This subroutine estimates H-P smoothered data.
Subroutine HP_TS

if !GapChoice = 1 then
	hpf(!NTS) GDP GDP_HP				' This is the trend;
	genr GDP_cyc = GDP - GDP_HP
	stomna(GDP_HP,GDP_HP_vec)
	stomna(GDP_cyc,GDP_cyc_vec)
	!lastgap = GDP_cyc_vec(!lasto)/GDP_HP_vec(!lasto)
else
	!lastgap=!lastgap
endif

Endsub		' End of subroutine HP_TS.
'****************************************************************************************************************************************



'****************************************************************************************************************************************
' This subroutine calculates the VAR regression.
Subroutine GetVAR

' VAR will be estimated with GDP growth rates and not GDP levels, so we first need to estimate growth.
genr growth = dlog(GDP)
genr l_REER = log(REER)


' This is just to estimate automatically the number of time series in %x_var. It changes the previous
' !nTS.
group tmp1 {%x_var}
stom(tmp1,mtmp1)
!nTS = @columns(mtmp1)

if !ch_lag = 0 then
	!Nlags = !setlag
else
' This procedure estimates the number of lags in the VAR that produces the minimum Akaike for the system.
!firsttmp = 0
vector(!maxlags) checkAkaike
for !ivar = 1 to !maxlags

	' Estimating the VAR with Nlags and a constant.
	var first_var.ls 1 !ivar {%x_var}  @

	!akaike=first_var.@aic				' Akaike Info. Criterion for the VAR system
	checkAkaike(!ivar) = !akaike		' This is just to check if Min Akaike is working

	if !firsttmp = 0 then
		!Nlags = !ivar
		!minak = !akaike
		!firsttmp = 1
	else 
		if !akaike < !minak then
			!Nlags = !ivar
			!minak = !akaike
		endif
	endif
	
next	' !ivar
endif
' Estimating the VAR with Nlags and a constant.
var first_var.ls 1 !Nlags {%x_var}  @


' Getting the covariance matrix of error terms.
sym res_var=first_var.@residcov	

' Estimating the Cholesky decomposition of the covariance matrix.
matrix cho_cov=@cholesky(res_var)

' Getting the coefficients matrix.
matrix coef=first_var.@coefmat	   ' need to get Bs separately

' Estimating the initial condition (VAR steady state).
' First I need to recover the B coefficients in an appropriate format, just as in the classic VAR specification..
for !i = 1 to !Nlags
	!temp = 0
	matrix(!nTS,!nTS) B{!i}
	for !j = 1 to !NTS
		for !k = 1 to !nTS
			B{!i}(!k,!j) = coef(!i+!temp,!k)
		next     ' !k
		!temp = !temp +!Nlags
	next     ' !j
next     ' !i

' Now I need to recover the constant.
vector(!nTS) C_var
for !ck1 = 1 to !nTS
	C_var(!ck1)= coef(!nTS*!Nlags + 1,!ck1)
next 	' !ck1

' Initial condition is then calculate, as defined in Ferrucci and Penalver (2003).
matrix(!nTS,!nTS) Bsum
for !ck2 = 1 to !Nlags
	Bsum = Bsum + B{!ck2}
next	' !ck2
matrix X_star = @inverse(@identity(!nTS) - Bsum)*C_var

'getting potential growth from X_Star
!potgrowth = (1 + X_star(!nTS-1))^4 - 1


Endsub		' End of subroutine GetVAR.
'****************************************************************************************************************************************



'****************************************************************************************************************************************
' This subroutine calculates the initial condition for the stochastic VAR.
Subroutine CalcInitialCond

' Two options are possible for initial conditions (others can be implemented in the future):
' (i) to use the steady-state condition as proposed in Ferrucci and Penalver (2003); or
' (ii) to use the last period for which there is data for all variables.

If !IniCond = 1 then
	' Initial condition using X_star. We need to create a vector of initial X-vars to be used in the forecasts. Initial 
	' condition must have the same number of lags as the number of forecasts.
	for !cc3 = 1 to !Nlags
		matrix X_lag{!cc3} = X_star
	next	' !cc3
else
	' Use the last observations.
	!nc = 0
	for !cc3 = !Nlags to 1 step -1
		if growth(!lasto - !nc) = NA or int_rate(!lasto - !nc) = NA then
			' In this case, there is missing observations in one of the lags. Then it will turn the steady-state case on
			' automatically, and shut the option 2 off.
			for !cc4 = 1 to !Nlags
				matrix X_lag{!cc4} = X_star
			next	' !cc4
			exitloop
		else
			vector(!nTS) X_lag{!cc3}
			X_lag{!cc3}(1) = growth(!lasto - !nc)
			X_lag{!cc3}(2) = int_rate(!lasto - !nc)
			!nc = !nc + 1
		endif
	next	' !cc3
endif

Endsub		' CalcInitialCond.
'****************************************************************************************************************************************



'****************************************************************************************************************************************
' This subroutine simulates stochastic forecasts for VAR.
Subroutine StochVAR

' Calculating forecasts using the autoregressive expression proposed in Ferrucci and Penalver (2003).
for !ck4 = 1 to !nSim
	matrix(!nTS,!nfor) X_for{!ck4}
	for !ck5 = 1 to !nfor
		matrix sumx = C_var
		for !ck6 = 1 to !Nlags
			sumx = sumx + B{!ck6}*X_lag{!ck6}
		next	' !ck6

		' Generating random errors from normal distributions with zero mean and same variance as in the VAR var-covar
		' matrix. Each vector of errors for each period is then pre-multiplied by the Cholesky factorization of VAR var-covar.
		matrix(!nTS) errs_var
		for !ck7 = 1 to !nTS
			series temp1 = nrnd
			stom(temp1,err_t)
			errs_var(!ck7) = err_t(1,1)
		next	' !ck7

		' This is the forecast.
		matrix mtemp = sumx + cho_cov*errs_var

		' Have to update the X_lags now for next forecasts.
		for !ck8 = !Nlags to 2 step -1
			!tt = !ck8 - 1
			X_lag{!ck8} = X_lag{!tt}
		next	' !ck8
		!ttt = 1
		X_lag{!ttt} = mtemp
		for !ck9 = 1 to !nTS
			X_for{!ck4}(!ck9,!ck5) = mtemp(!ck9)
		next	' !ck9

	next	' !ck5

next	' !ck4

' This procedure generates forecast quarters and years that are consistent with the last quarter of data (for all TS).
matrix(!nfor,1) obsfor
matrix(!nfor,1) Qtfor
matrix(!nfor,1) Yearfor
if !Lastq < 4 then
	Qtfor(1) = !Lastq + 1
	Yearfor(1) = !Lasty
else
	Qtfor(1) = 1
	Yearfor(1) = !Lasty + 1
endif

for !ck10 = 2 to !nfor
	if Qtfor(!ck10 - 1) < 4 then
		Qtfor(!ck10) = Qtfor(!ck10 - 1) + 1
		Yearfor(!ck10) = Yearfor(!ck10 - 1)
	else
		Qtfor(!ck10) = 1
		Yearfor(!ck10) = Yearfor(!ck10 - 1) + 1
	endif
next	' !ck10

obsfor(1) = !lasto + 1
for !ck11 = 2 to !nfor
	obsfor(!ck11) = obsfor(!ck11 - 1) + 1
next	' !ck11

Endsub		' End of subroutine StochVAR.
'****************************************************************************************************************************************



'****************************************************************************************************************************************
' This subroutine gets annual forecasted X_var data for forecasting stock of debt.
Subroutine Annualize

' We have to be careful here. For forecasting debt, we will use the first expression in Ferrucci and Penalver (2003), in
' which current total debt level depends on lag total debt, real GDP growth rate, real interest rate, and primary surplus.
' There is another version for debt forecast, which disentangles debt into foreign and domestic debt and which includes
' other variables such as interest rate in foreign and domestic currency and currency depreciation.
' For the simpler version, annualized GDP growth rates have to be compounded for annualization. The other variables
' will be used as levels and we just need to get the 4th quarter data.
' For the second version, currency depreciation would have to be comounded as well, but the same rationale used for
' GDP growth rates can be easily extended for currency depreciation.

' This matrix stores original data used in the VAR estimation. Some historical data is used, when the last historical
' period is not quarter 4.
group x_vard {%x_var}
stomna(x_vard,x_varm)

for !cn1 = 1 to !nSim
	!cc = 1
	matrix(!nTS,!Nfor_a) X_for_a{!cn1}
	for !cn2 = 1 to !nfor

		if Qtfor(!cn2) = 4 then
			for !kc = 1 to !nTS

				' The following procedure works fine if there are no missing values (NA) for growth rates prior to the forecast
				' period.
				if annualize(!kc) = 1 then	' Will compound it

					if obsfor(!cn2) - !lasto <= 4 then

						if !lastq = 4 then
							X_for_a{!cn1}(!kc,!cc) = (1 + X_for{!cn1}(!kc,!cn2 - 3))*(1 + X_for{!cn1}(!kc,!cn2 - 2))*(1 + X_for{!cn1}(!kc,!cn2 - 1))*(1 + X_for{!cn1}(!kc,!cn2)) - 1
						else 
							if !lastq = 3 then
								X_for_a{!cn1}(!kc,!cc) = (1 + X_for{!cn1}(!kc,!cn2 - 2))*(1 + X_for{!cn1}(!kc,!cn2 - 1))*(1 + X_for{!cn1}(!kc,!cn2))*(1 + x_varm(!lasto,!kc)) - 1
							else 
								if !lastq = 2 then
									X_for_a{!cn1}(!kc,!cc) = (1 + X_for{!cn1}(!kc,!cn2 - 1))*(1 + X_for{!cn1}(!kc,!cn2))*(1 + x_varm(!lasto,!kc))*(1 + x_varm(!lasto-1,!kc)) - 1
								else 
									if !lastq = 1 then
										X_for_a{!cn1}(!kc,!cc) = (1 + X_for{!cn1}(!kc,!cn2))*(1 + x_varm(!lasto,!kc))*(1 + x_varm(!lasto-1,!kc))*(1 + x_varm(!lasto-2,!kc)) - 1
									endif
								endif
							endif
						endif
						
					else
						X_for_a{!cn1}(!kc,!cc) = (1 + X_for{!cn1}(!kc,!cn2 - 3))*(1 + X_for{!cn1}(!kc,!cn2 - 2))*(1 + X_for{!cn1}(!kc,!cn2 - 1))*(1 + X_for{!cn1}(!kc,!cn2)) - 1
					endif

				else		' Will average it

					' Have to find the denominator, depending upon whether there any missing observations in the data.
					' We know the last historical period prior to forecast have data for all variables (this is the way this
					' program has been set up. We just need to confirm if the previous two have missing values or not
					if x_varm(!lasto-1,!kc) = NA then
						!den1 = 3
						if x_varm(!lasto-2,!kc) = NA then
							!den2 = 2
						endif
					else 	' There is data for x_varm(!lasto-1,!kc)
						if x_varm(!lasto-2,!kc) = NA then
							!den2 = 3
						else
							!den1 = 4
							!den2 = 4
						endif
					endif

					if obsfor(!cn2) - !lasto <= 4 then
						if !lastq = 4 then
							X_for_a{!cn1}(!kc,!cc) = ( X_for{!cn1}(!kc,!cn2 - 3) + X_for{!cn1}(!kc,!cn2 - 2) + X_for{!cn1}(!kc,!cn2 - 1) + X_for{!cn1}(!kc,!cn2))/4
						else 
							if !lastq = 3 then
								X_for_a{!cn1}(!kc,!cc) = ( X_for{!cn1}(!kc,!cn2 - 2) + X_for{!cn1}(!kc,!cn2 - 1) + X_for{!cn1}(!kc,!cn2) + x_varm(!lasto,!kc))/4
							else 
								if !lastq = 2 then
									X_for_a{!cn1}(!kc,!cc) = ( X_for{!cn1}(!kc,!cn2 - 1) + X_for{!cn1}(!kc,!cn2) + x_varm(!lasto,!kc) + x_varm(!lasto-1,!kc))/!den1
								else 
									if !lastq = 1 then
										X_for_a{!cn1}(!kc,!cc) = ( X_for{!cn1}(!kc,!cn2) + x_varm(!lasto,!kc) + x_varm(!lasto-1,!kc) + x_varm(!lasto-2,!kc))/!den2
									endif
								endif
							endif
						endif
					else
						X_for_a{!cn1}(!kc,!cc) = ( X_for{!cn1}(!kc,!cn2 - 3) + X_for{!cn1}(!kc,!cn2 - 2) + X_for{!cn1}(!kc,!cn2 - 1) + X_for{!cn1}(!kc,!cn2))/4
					endif


				endif

			next	' !kc

			!cc = !cc + 1

		endif


	next	' !cn2

next	' !cn1

Endsub		' End of subroutine Annualize.
'****************************************************************************************************************************************



'****************************************************************************************************************************************
' This subroutine forecasts output gap, primary surplus, and debt.
Subroutine ForGapPSDebt

' This subroutine has been implemented assuming the linear model for PS has been previously estimated and
' parameters have been set in the Setup subroutine above, along with all other parameters. It is possible to estimate
' the PS regression here and use the estimated coefficients for forecasting PS. Ideally this regression would be
' estimated every time a new debt forecast is calculated.


if !PS_Choice = 0 then
	vector(!nfor_a) tempz
	PS_Shocks = tempz
endif

matrix(1,1) tttemp
tttemp(1,1) = !Nfor_a
for !cgap1 = 1 to !nSim
	!old_gdp = gdp(!lasto)
	!old_gap = !lastgap
	!old_debt = !Curdebt
	vector(!Nfor_a) forgap{!cgap1}
	vector(!Nfor_a) forps{!cgap1}
	vector(!Nfor_a) fordebt{!cgap1}
	for !cgap2 = 1 to !Nfor_a
		if x_for_a{!cgap1}(1,!cgap2) <> 0 then
				
			forgap{!cgap1}(!cgap2) = (1 + !old_gap)*(1+x_for_a{!cgap1}(!nTS-1,!cgap2))/(1+!potgrowth) - 1	' New forecasted output gap

			'This part is just to estimate the random error term
			series temp1 = nrnd
			stom(temp1,err_ps)
			series temp2 = !var_ps
			stom(temp2,ps_var)
			!errps_var = err_ps(1,1)*@sqrt(ps_var(1,1))

			forps{!cgap1}(!cgap2) = !c_ps + PS_Shocks(!cgap2) + !b1_ps*forgap{!cgap1}(!cgap2) + !b2_ps*!old_debt + !errps_var

			if !debt_choice = 0 then
				fordebt{!cgap1}(!cgap2) = (1 + X_for_a{!cgap1}(!nTS-1,!cgap2) - X_for_a{!cgap1}(!nTS-2,!cgap2))*!old_debt - forps{!cgap1}(!cgap2)

			else
				' New debt expression

				' First: foreign component
				if !cgap2 = 1 then
					!app_REER=X_for_a{!cgap1}(!nTS,!cgap2) - x_varm(!lasto,!nTS)
				else
					!app_REER=X_for_a{!cgap1}(!nTS,!cgap2) - X_for_a{!cgap1}(!nTS,!cgap2-1)
				endif

				!for_comp = (((1 - !app_REER)*(1 + X_for_a{!cgap1}(!nTS-3,!cgap2)))/(1 + X_for_a{!cgap1}(!nTS-1,!cgap2)))*debt_s(!cgap2)*!old_debt
	
				' Second: domestic component
				!dom_comp = ((1 + X_for_a{!cgap1}(!nTS-2,!cgap2))/(1 + X_for_a{!cgap1}(!nTS-1,!cgap2)))*(1-debt_s(!cgap2))*!old_debt

				' Bringing them together
				fordebt{!cgap1}(!cgap2) = !for_comp + !dom_comp - forps{!cgap1}(!cgap2)

			endif

			!old_gap = forgap{!cgap1}(!cgap2)					' Updated output gap for next period
			!old_debt = fordebt{!cgap1}(!cgap2)					' Updated stock of debt for next period
		endif
	next	' !cgap2
next	' !cgap1

if !outp = 1 then
	!ini1 = !lasty + 1
	!fin1 = !lasty + !Nfor_a
	range !ini1 !fin1
	smpl !ini1 !fin1
	for !cgap3 = 1 to !nSim
		mtos(fordebt{!cgap3},TSForecastDebt{!cgap3})
		if !cgap3 = 1 then
			group fordebt TSForecastDebt{!cgap3}
		else
			fordebt.add TSForecastDebt{!cgap3}
		endif
	next	' !cgap3

	write(t) {%OutputPath}\debt_forecast.txt fordebt

endif


Endsub		' End of subroutine ForGapPSDebt.
'****************************************************************************************************************************************



'****************************************************************************************************************************************
' This subroutine estimates trends for current and for forecasted debt.
Subroutine GetTrend

' This subroutine has been implemented as a first step for generating the fan charts (Bank of England).
' First, we need to generate the trend for each forecasted debt path. Then, we estimate lines.

matrix(!nSim,2) B1_fordebt

for !ktrend1 = 1 to !nSim

	' The series for estimating the trend will include current debt level (as a fraction of GDP) plus all
	' forecasted debt
	stom(TSForecastDebt{!ktrend1},fordtemp{!ktrend1})
	vector(!Nfor_a + 1)  nfordtemp{!ktrend1}
	vector(!Nfor_a + 1) obsT
	for !ktrend0 = 1 to !Nfor_a + 1
		if !ktrend0 = 1 then
			nfordtemp{!ktrend1}(!ktrend0) = !curdebt
			obsT(!ktrend0) = 1
		else
			nfordtemp{!ktrend1}(!ktrend0) = fordtemp{!ktrend1}(!ktrend0-1)
			obsT(!ktrend0) = obsT(!ktrend0 - 1) + 1
		endif
	next	' !ktrend0
	range !lasty !fin1
	smpl !lasty !fin1
	mtos(nfordtemp{!ktrend1},nTSForecastDebt{!ktrend1})
	mtos(obsT,nobsT)

	' Estimate the trend using a simple OLS: for_debt(t) = C + B1*t, t =1, ..., T
	equation eq{!ktrend1}.ls nTSForecastDebt{!ktrend1} c nobsT
	vector coef_trend=eq{!ktrend1}.@coefs 	' Getting the coefficients matrix.
	B1_fordebt(!ktrend1,2) = !ktrend1
	B1_fordebt(!ktrend1,1) = coef_trend(2)

	' Now generates the TS with the estimated coefficients and starting at the initial debt
	vector(!Nfor_a+1) debt_trend{!ktrend1}
	debt_trend{!ktrend1}(1) = !curdebt
	for !ktrend2 = 2 to !Nfor_a + 1
		debt_trend{!ktrend1}(!ktrend2) = debt_trend{!ktrend1}(!ktrend2-1)+ B1_fordebt(!ktrend1,1)
	next	' !ktrend2

next	' !ktrend1

' Sorting the debt forecasts for estimating the deciles to be plotted in the fan chart
!fin2 = !lasty + !nSim - 1
range !lasty !fin2
smpl !lasty !fin2
mtos(B1_fordebt,B1_fordebtg)
sort B1_fordebtg
stom(B1_fordebtg,SortedB1_fordebtg)


' Getting the deciles
for !kdec1 = 10 to 90 step 10

	vector dec_position{!kdec1} = SortedB1_fordebtg(@round(!kdec1*!nSim/100),2)

next	' !kdec1


range !lasty !fin1
smpl !lasty !fin1

for !kdec2 = 10 to 90 step 10
	!temp = dec_position{!kdec2}(1)
	mtos(debt_trend{!temp},dec_trend{!kdec2})

	' Putting all deciles in one single group
	if !kdec2 = 10 then
		group gtrenddebt_dec dec_trend{!kdec2}
	else
		gtrenddebt_dec.add dec_trend{!kdec2}
	endif
next	' !kdec2

write(t) {%OutputPath}\debt_trend.txt gtrenddebt_dec


range !debtfirstyr !fin1
smpl !debtfirstyr !fin1

Read(t=txt, d=t) {%XSourceData2} hist_debt					' Reads historical debt data

Graph Fan_Chart.line hist_debt gtrenddebt_dec
Fan_Chart.addtext(t) Debt Forecasts: {%Country_Name}					' Graph title
Fan_Chart.addtext(2.5,-.3,l) Debt (% of GDP)						' Y-axis label
Fan_Chart.addtext(1.8,3.5) Years										' X-axis label

Fan_Chart.options linepat												' Set chart options
Fan_Chart.legend position(r)											' Set legend position

' Following lines format the series, starting by historical debt and going through the deciles forecasts.
Fan_Chart.setelem(1) lpat(solid) lcolor(blue)					
Fan_Chart.setelem(2) lpat(dash1) lcolor(blue)	
Fan_Chart.setelem(3) lpat(dash7) lcolor(red)	
Fan_Chart.setelem(4) lpat(solid) lwidth(2) lcolor(green)	
Fan_Chart.setelem(5) lpat(solid) lwidth(3) lcolor(blue)	
Fan_Chart.setelem(6) lpat(solid) lwidth(1) lcolor(red)	
Fan_Chart.setelem(7) lpat(solid) lwidth(3) lcolor(blue)	
Fan_Chart.setelem(8) lpat(solid) lwidth(2) lcolor(green)	
Fan_Chart.setelem(9) lpat(dash7) lcolor(red)	
Fan_Chart.setelem(10) lpat(dash1) lcolor(blue)	

Endsub		' End of subroutine GetTrend.
'****************************************************************************************************************************************



'****************************************************************************************************************************************
' This subroutine estimates the probability that the debt may fall above or below certain thresholds, pre-defined in the SetUp
' Subroutine by user.
Subroutine GetProb


vector(!Nfor_a) probup
vector(!Nfor_a) probdown

for !cprob1 = 1 to !Nfor_a

	for !cprob2 = 1 to !nSim

		if fordebt{!cprob2}(!cprob1) > !upperdebt then
			probup(!cprob1) = probup(!cprob1) + 1
		endif

		if fordebt{!cprob2}(!cprob1) < !lowerdebt then
			probdown(!cprob1) = probdown(!cprob1) + 1
		endif

	next	'!cprob2

next	'!cprob1

vector fprobup = probup/!nSim
vector fprobdown = probdown/!nSim

if !outp = 1 then
	mtos(fprobup,TSfprobup)
	mtos(fprobdown,TSfprobdown)
	write(t) {%OutputPath}\Upper_thres_probab.txt TSfprobup
	write(t) {%OutputPath}\Lower_thres_probab.txt TSfprobdown
endif

Endsub		' End of subroutine GetProb.
'****************************************************************************************************************************************



'****************************************************************************************************************************************
' This subroutine estimates percentage cutoffs on the forecasted debt, for each period.
Subroutine CalcStat

'range 1 !Lastobs
'smpl 1 !Lastobs
' First we create vectors to store all simulated values for debt for each forecast period.
for !cstat1 = 1 to !nSim
	!ctemp = 0
	for !cstat2 = 1 to !Nfor_a
		if x_for_a{!cstat1}(1,!cstat2) <> 0 then
			vector(!nSim) DebtStat{!cstat2}
			DebtStat{!cstat2}(!cstat1) = fordebt{!cstat1}(!cstat2)
			!ctemp = !ctemp + 1
		endif
	next	' !cstat1
next	' !cstat2
'stop
' Now we sort (ascending) each of the DebtStat vectors, so as to be able to pick any percentage cut-off.
!tmptest=(1990+!nSim-1)
range 1990 !tmptest
smpl 1990 !tmptest
for !cstat3 = 1 to !ctemp
	vector(!nSim) SortedDebt{!cstat3}
	mtos(DebtStat{!cstat3},forecast_debt_period{!cstat3})
	forecast_debt_period{!cstat3}.hist
	sort(d) forecast_debt_period{!cstat3}
	stom(forecast_debt_period{!cstat3},SortedDebt{!cstat3})
next	' !cstat3

vector(!Nfor_a) percDebt
for !cstat4 = 1 to !ctemp
	' !perc of the time the debt is below percdebt variable.
	percdebt(!cstat4) = SortedDebt{!cstat4}(@ceiling((1-!perc)*!nSim))
next	' !cstat4

if !outp = 1 then
	mtos(percdebt,TSpercdebt)
	write(t) {%OutputPath}\debt_percentile.txt TSpercdebt
endif

Endsub		' End of subroutine CalcStat.
'****************************************************************************************************************************************

