'*******************************************************************************************************
' Evaluation of Exchange Rate Forecasts 
' Program: exchange6.prg
' Workfile: nexra2.wf1, nexra4.wf1
' First version: 10.29.08
' Last update: 07.16.10
' Roberto Duncan
'*******************************************************************************************************
'
' General Comments:
'  This program is designed to work in Eviews 6 by using workfiles nexra2.wf1 or nexra4.wf1. 
'  Use nexra2 to obtain results for forecast evaluation. In particular, 
' CW, UTheils stats, p-values and max t-stats (as in, e.g., emw23 or emw26.xls). 
' Also, it can be used to forecast the Euro using four models (see emw20.pdf)
' and data up to 2007.IV (e.g., as in emw28.xls).
' Use nexra4 to forecast the Euro using the most updated exchange rate series and model 1, 
' the model with factors only (e.g. as reported in "all models" sheet of emw norges2.xls)
' Both stats and forecasts are obtained for several forecast horizons. Note that we make Euro 
' forecasts by averaging the 8 individual forecasts. Main output matrices are "a_ucw_s" and 
' "a_feuro", also saved as Excel sheets (see section 7).
'
' Workfiles:
' -Nexra 2 contains 4 series for each country: m_i, p_i, s_i, and y_i, for i=1,2,...,18. The letters stand
' for money supplies, price  indexes, nominal exchange rates (foreig currency per US$), and real
' outputs, respectively; all expressed in logs for the 1973.I-2007.IV period (140obs). The  countries   
' are: AUS, CAN, DEN, GBR, JPN, KOR, NOR, SWE, SWI, AUT, BEL, FRA, GER, ESP, ITA, 
' FIN, NET, and USA. Note that s_18=0 for every period since i=18 corresponds to USA. Aside from 
' this, note also that Euro area countries are placed towards the end of the list.
' -Nexra 4 contains the same data (same series. ordering, etc.) but exchange rates series are
' updated up to 2010.II (150obs). Note: last quarter is as of May 2010.
'
' Related Documents/Files: 
' Results of forecast evaluation: emw23.xls, emw25.xls, emw28.xls, among others
' Results of Euro forecasts: emw norges2.xls, etc.
' Requests, background, models, etc: “inf3.pdf, emw20.pdf, emw21.pdf, etc
'
'*******************************************************************************************************
'
' Contents: 
' The code has 7 sections:
' 1. Subroutines
' 	Four subroutines (HP filtering, long run variance, generation of random draws,  
'  	and p-values for one-sided test), called in sections 3 and 6, are constructed.
' 2. Main Parameters
' 	Sets the main parameters of the simulations and forecasts (basically type of 
' 	sample and model, number of factors, whether Euro forecast is desired or not, etc.)
' 3. Constructing Matrices and Series
'	Constructs output matrices and series used in the the estimation of models
' 	 and Euro forecasts, performs seasonal adjustment and starts main loop for each
' 	forecast horizon.
' 4. Recursive Estimation and Forecasts
' 	Starts the recursive estimation; creates factor group and estimates factor loadings
' 	 (using unrotated solution), estimates one of four types of models (only factors, 
'	with Taylor Rule, monetary fundamentals, or PPP).
' 5. Euro Area Forecast
'	Euro forecasts for several horizons using simplest factor model and full sample
' 6. Forecast Evaluation
' 	Constructs prediction errors, Theil's U stats, Clark-West stats for univariate
' 	and multivariate cases, max t-stat and p-values
' 7. Excel Output
'	Saves output in an Excel spreadsheet in the specified directory.
'
'*******************************************************************************************************
'
' Instructions:
' 1- Open workfile (nexra2.wf1 or nexra4.wf1) by clicking on "File" and then "Eviews Workfile".
' 2- Choose parameter values in section 2 below. Set
' 	%S to "F", "E" or "L" to control if full, early (pre Euro) or late (post Euro) sample is used
' 	M to 1, 2, 3, or 4 to control whether "only factors" (F-s), "Taylor Rule" (F-s,tr), "monetary" 
'		(F-s, mf), or "PPP" (F-s, ppp) model is used.
' 	F to fix the number of factors (F=1,...,4) in the model chosen
' 	%Euro to "Yes" if forecast of Euro is also required. "No" does not generate Euro forecasts.
'		Notes: it works only with %S="F". If nexra4.wf1 is used then set M=1 (model with 
'		only factors) unless series for prices, ouput, money supplies are updated. 
'
'	ND to fix the number of random draws for critical values and max t-stat. Note: ND=10000
'		was used to obtain results reported in emw23.xls.
'
'	Note: There are some parameters that do not need to be changed unless a major
'		change to the code has to be done or a particular result is needed:
'	N1 is fixed equal to the number of currencies or exchange rates (currently 17)
'	R is already set equal to the number of observation in range of the workfile. For 
'		example, in nexra2.wf1, R will be 140.
'	%fp is set to "1973q1" for the first period in the workfile range.
'	mvsa is set to 3, the # of missing values due to seasonal adjustment
'	T0 is 55, the place of the first forecast period (e.g.73q1+t0=86.1)
'	PFE=R, the # of quarters since 73q1 to use a forecast base of Euro
' 	The vector "horizon" that determines the forecast horizons (currently: 1,4,8,12,16)
' 3- Change the location and name of directory in section 7 at the end of this code in order
'	to let the code save the output  as needed.
' 4- Press [Run]. The program takes a few seconds to run and store/save results.
'
' Main results related to forecast evaluation (CW, UTheils, stats; max-stats) can be found either 
' in the vector named "a_ucw_s" in the workfile or in the Excel sheet saved in the directory
' specified in previous step. The file name will be "stats_M_F.xls" for M and F set in step 2.
' If %Euro="Yes" was chosen in step 2, the forecasts of the Euro are stored as a vector of the
' workfile with the name "a_feuro".
'
'*******************************************************************************************************
'*******************************************************************************************************
'
' 1. Subroutines:
' S1. Recursive HP filter (rhpfilter; called in section 3)
' S2. Long-run variance (gam0; called in section 6)
' S3. Random draws from a N(0,Omega) (rdraw; called in section 6)
' S4. P-value for one-sided test based on simulated critical values (pvals; called in section 6)
'
'********************************************************************************************************
'
' S1. Recursive HP filter
' Objective: Detrends series recursively using the Hodrick-Prescott filter
' Note: It is a slightly modified version of Nelson Mark's subroutine
' Inputs:		bigt  scalar, last obs # of the recursion 
'	    	x   raw/original series
'    		t1  scalar, starting # of obs for the recursion (minimum 4)
'  Output:'    	xf  detrended or filtered series.

 subroutine local rhpfilter(scalar bigt, series x, scalar t1, series xf)
 poff                               	' no print 
 %1 = @otod(t1)                    ' convert obs no to date
 smpl @first %1                    ' start HP with the first t1 obs
 do x.hpf hptr                  	 ' HP (uses default smoothing parameter)
 xf = x - hptr                     	' detrended series
 for !j= t1+1 to bigt              	' recursion
   %1 = @otod(!j)                  
   smpl @first %1                  
   do x.hpf hptr                       
   xf(!j)=x(!j)-hptr(!j)            
 next
 smpl @all                        	 ' restore data range
 endsub
'
'********************************************************************************************************
'
' S2. Long-run variance 
' Objective: Estimates long run variance matrix for a set of variables (e.g. forecasts)
' Inputs:		P   scalar, maximum number of forecasts
'		k   scalar, forecast horizon
'		N1 scalar, number of countries or individuals
'		Dbig (N1xP-k+1) matrix with P-k+1 forecasts of N1 series
' Output:		Gama0 N1xN1 autocovariance matrix of order j=0

  subroutine local gam0(scalar P, scalar k, scalar N1, matrix Dbig, matrix Gama0)
  poff                          ' no print 
  matrix (N1,N1) Gama0
  for !j = 1 to P-k+1
	vector v_{!j} = @columnextract(dbig,!j)
	matrix m_{!j}=v_{!j}*@transpose(v_{!j})
	matrix Gama0=Gama0+m_{!j}
  next
  matrix Gama0=(1/(P-k+1))*Gama0
  endsub
'
'********************************************************************************************************
'
' S3. Random Draws from a N(0,Omega)
' Objective: Generates random draws from a Normal distribution with correlation matrix Omega
' Comment: Uses a fixed random seed
'  Inputs:		ND   scalar, number of draws
'		NE  scalar, number of elements
'		ome matrix, correlation matrix
'  Output: 	z matrix, NDxNE matrix with random draws
'		rankmaxz vector, vector with ND ranked maximum values of the NE random draws

 subroutine local rdraw(scalar nd, scalar ne, matrix ome, matrix z, vector rankmaxz)
 poff                               ' no print 
 matrix(nd,ne) un
 rndseed 123456
 for !j=1 to nd
  for !h=1 to ne
   un(!j,!h)=@rnorm
  next
 next
 sym omex = ome
 matrix sqome = @cholesky(omex)
 matrix z=@transpose(sqome)*@transpose(un)
 vector maxz = @cmax(z)
 vector rankmaxz = @sort(maxz,"a","i")
 endsub
'
'********************************************************************************************************
'
'  S4. P-value for one-tailed test based on simulated critical values
'  Inputs:		nd   scalar, number of elements in vector of simulated critical values
'		zhat  vector, 1x1 with t-statistic
'		rz vector, vector with nd critical values 
'  Output: 	p_val 1x1 vector (fraction of the nd values in rz that are larger than zhat

 subroutine local pvals(scalar nd, vector zhat, vector rz, vector p_val)
 poff                               ' no print 
 vector (nd) auxz
 scalar zhataux=zhat(1)
 for !j=1 to nd
  auxz(!j)=rz(!j)>zhataux
 next 
 vector sumauxz=@csum(auxz)
 vector p_val=sumauxz/nd
 endsub
'
'********************************************************************************************************
'********************************************************************************************************
'
' 2. Main Parameters
'
'********************************************************************************************************
' Set parameter values for %S, M, F, %EUR, ND
'
%S="F"				' E=early sample (pre Euro), L=late smpl (post), F=full smpl
scalar M=1			' model number  1,...,4, or (F-s), (F-s,tr), (F-s,mf), (F-s,ppp)
scalar F=2	     		' number of factors (F=1,...,4)
%EUR="Yes"			' "Yes" if forecast of Euro is needed; it works only with %S=F
scalar ND=10000 		' number of random draws for critical values and max t stat
'
'********************************************************************************************************
' Parameters N1, R, %fp, mvsa, T0, PFE or forecast horizons can be changed if a major 
' modification to the code or a particular result is needed. 
'
scalar N1 = 17             		' number of currencies or exchange rates except USA
scalar R=@obsrange		' number of obs in workfile range (e.g. nexra2.wf1: R=140)
%fp="1973q1"			' first period in workfile range
scalar mvsa=3			' number of missing values due to seasonal adjustment
scalar T0 = 55             		' first forecast period (e.g. 73q1+t0=86.1)
scalar PFE=R			' number of quarters since 73q1 to use forecast base of Euro

vector(5) horizon			' number of forecast horizons
horizon(1)=1			' forecast horizons; currently/usually 1, 4, 8, 12, 16 quarters
horizon(2)=4
horizon(3)=8
horizon(4)=12
horizon(5)=16
scalar rh=@rows(horizon)
'
'********************************************************************************************************
' To determine the maximum number of forecasts for the sample chosen
' Hard coded numbers: P=49 for early sample (pre Euro); =35 for late sample (post Euro).
' NDE=No obs to Drop at End of sample (e.g. nde=39 to start 98q1)
' P is max No of forecasts when forecast horizon is 1, equals 49 and 35 for early and late sample
'
if %S="F" then
scalar NDE=0
scalar P=R-t0-1-NDE     
endif

if %S="E" then
scalar NDE=32
scalar P=49	     		 
endif

if %S="L" then
scalar NDE=0
scalar P=35	     		 
endif
'
'********************************************************************************************************
'********************************************************************************************************
'
' 3. Constructing matrices and series
'
'********************************************************************************************************
' Matrices that store main results (a_ucw_s and a_feuro),
' This section also starts main loop for each forecast horizon, adjusts series (output, prices, money)
' for seasonal components, filters output series creating, 
' Notation:
' "a_ucw_s" stores UTheil, Clark-West stats, p-values and max t-stats
' Its first N1 rows for the U and CW stats of N1 countries, the N1+1th row contains only zeros, 
' the last two rows contain p-values and max t-stats for each forecast horizon
' See also last section for a description of the contents of a_ucw_s.
' "a_feuro" stores forecast of Euro for each forecast horizon
' ysa, psa, and msa denote seasonal adjusted series for output, prices, money balances
' pi denotes inflation rates; gap is detrende ysa
' nexra is group of N1 exchange rates; series shat(e,euro) are used to store forecasts
' nf = (auxiliary) vector with factor numbering
' Hard coded number in call rhpfilter: t1=7= 3 obs lost (seasonal adj.)+4 (min # obs to HP filtering) 
'  
matrix(N1+3,2*rh) a_ucw_s

if  %EUR="Yes" then		' only if Euro forecats is required
vector(rh) a_feuro
endif

for !hh = 1 to rh 		'main loop for each forecast horizon starts here
 scalar k=horizon(!hh)
 scalar she=0

	smpl %fp+mvsa  @last       
	for !j = 1 to N1+1
 		'seasonal adjustment of series (output, prices, money)
		' by averaging the last 4 quarters
  		series ysa_!j =(y_!j+y_!j(-1)+y_!j(-2)+y_!j(-3))/4
  		series psa_!j =(p_!j+p_!j(-1)+p_!j(-2)+p_!j(-3))/4 
  		series msa_!j =(m_!j+m_!j(-1)+m_!j(-2)+m_!j(-3))/4
  		'seasonal adjustment of inflation rates
   		series pi_!j = psa_!j-psa_!j(-1)
  		'generating output gap for model with Taylor rule
  		series x = ysa_!j		'original/raw series                  
  		series gap_!j = 0		'detrended series
  		scalar bigt = @obs(x)+mvsa 	'last obs # of the recursion  
  		call rhpfilter(bigt,x,7,gap_!j) '3 obs lost + 4 (min # obs to HP filtering)
	next

 smpl %fp @last   

 'forming group of exchange rates
 group nexra s_1
 for !j = 2 to N1  
  nexra.add s_{!j}
 next 
 series shate_{!hh}=0

 smpl 1999q1 @last-k
 series shateuro_{!hh}=0

 ' To store forecasts for N1 countries
 smpl @first @first   

 for !j = 1 to N1
      series shat_{!j} = 0  
 next

' Vector with factor numbering
  vector(F) nf
  for !f=1 to F
  nf(!f)=!f
  next

'********************************************************************************************************
'********************************************************************************************************
'
' 4. Recursive Estimation and forecasts
'
'********************************************************************************************************
' This section starts the recursive estimation; creates factor group and estimates and
' stores factor loadings (using unrotated solution), constructs F-s type regressors, estimates
' one of the four models to be used for forecasting
' Method for computing the factor score coefficient matrix: Green (coef=“green”). For furthe detail
' see Eviews 6 Command Reference, page 115.
' Notation:
' te= # obs to the upper limit of estimation loop for each sample
' mv=# missing values in lhs variable (for inflation work only)
' FS1,..., FS4 are factor series names
' "loads" stores factor loadings 
' Hard coded numbers in c(600+!j) for !j=1,..,N1, c(650), c(650) are the place of some coefficient 
' estimates stored transitorily in the vector "c' of the workfile.

 smpl @all
 scalar mv=R-@obs(pi_1)

 if %S="E" then
 scalar te=t0+P-1				'
 endif
 
 if %S="F" then
 scalar te=R-2-(k-1) 
 endif

 if %S="F" and %EUR="Yes" then
 scalar te=R-2-(k-1)+1 
 endif

 if %S="L" then
 scalar t0=104
 scalar te=t0+P-k
 endif

'********************************************************************************************************
'   Estimating factor series and loadings, saving estimated factor score series in the workfile

 for !t =  t0 to te 		' loop starts up to upper limit te

   smpl %fp+mv  %fp+!t   ' recursive updating
 
   factor factor02.ml(n=F) nexra 'creating factor group (n factors)
  
  ' naming factors (no more than F factors, see "Main parameters" above)
   if F=1 then
   factor02.factnames FS1      
   endif

   if F=2 then
   factor02.factnames FS1  FS2    
   endif

   if F=3 then
   factor02.factnames FS1  FS2  FS3
   endif

   if F=4 then
   factor02.factnames FS1  FS2  FS3  FS4
   endif
  
  factor02.makescores(coef=green,n=factors_outs) nf  'extracting factor series
  'note that the method for computing factor score coefficient matrix is Green's 

  matrix loads = factor02.@loadings   'stores factor loadings (using unrotated solution)

'********************************************************************************************************
'   Creating a pool, constructing regressors F(it)-s(it) for 1,...,F factors, i=1,...,N1 
'   currencies/countries. It formulates ans estimates one of the four models per value 
'   of M in Main Parameters

   system POOL6          

   for !j = 1 to N1  
   
   series Fsx_{!j} = - s_{!j} 
     for !f=1 to F
	series Fsx_{!j} =Fsx_{!j}+ loads(!j,!f)*FS{!f}  'Regressors F(it)-s(it) 
     next

    series Fs_{!j} =Fsx_{!j}

     	if M=1 then
    	POOL6.append  s_!j-s_!j(-k) = c(600+!j) + c(650)*Fs_{!j}(-k)
    	endif
     	if M=2 then
	series xtr_{!j} =1.5 * (pi_{!j} - pi_18)+0.5*(gap_{!j} -gap_18)
    	POOL6.append  s_!j-s_!j(-k) = c(600+!j) + c(650)*Fs_{!j}(-k)+c(660)*xtr_{!j}(-k)
    	endif
     	if M=3 then
	series xmm_{!j} = msa_{!j} - msa_18 - ( ysa_{!j} - ysa_18 ) - s_{!j}
    	POOL6.append  s_!j-s_!j(-k) = c(600+!j) + c(650)*Fs_{!j}(-k)+c(660)*xmm_{!j}(-k)
    	endif	 
     	if M=4 then
	series xppp_{!j} = psa_{!j} - psa_18  - s_{!j}
    	POOL6.append  s_!j-s_!j(-k) = c(600+!j) + c(650)*Fs_{!j}(-k)+c(660)*xppp_{!j}(-k)
    	endif	 
   next

  do POOL6.ls

'********************************************************************************************************
'  Storing last obs in series and generating forecasts

  smpl %fp + !t  %fp + !t    
   for !j = 1 to N1      	 
	if M=1 then
    	series sto{!j} = c(600+!j) + c(650)*Fs_{!j}
	endif
	if M=2 then
    	series sto{!j} = c(600+!j) + c(650)*Fs_{!j}+c(660)*xtr_{!j}
	endif
	if M=3 then
    	series sto{!j} = c(600+!j) + c(650)*Fs_{!j}+c(660)*xmm_{!j}
	endif
	if M=4 then
    	series sto{!j} = c(600+!j) + c(650)*Fs_{!j}+c(660)*xppp_{!j}
	endif
    shat_{!j}=sto{!j}
    series shat_{!j}_{!hh}=shat_{!j}
   next
   
 next			' loop ends here

'********************************************************************************************************
'********************************************************************************************************
'
' 5. Euro area forecast
'
'********************************************************************************************************
' The forecasts are stored in "a_feuro", a vector in the workfile. Probably the second object stored.
' Check that %EUR="Yes" and %S="F" (it works mainly with full sample).
' If nexra4.wf1 is used then set M=1 (model with 
' only factors) unless series for prices, ouput, money supplies are updated in a new workfile. 
' Note that "/8" is because we make euro forecasts by averaging the 8 individual forecasts. 

   'Only for late sample (%S="L"). This is just for forecast as average of Euro area forecasts
   if %S="L" then
    	smpl 1999q1 @last-k
    	for !e=10 to N1
      		shate_{!hh}=shate_{!hh}+shat_{!e} 
    	next
    	series shateuro_{!hh}=shate_{!hh}/8	'average of Euro area forecasts
   endif

'********************************************************************************************************
   'Only for full sample (%S="F") and forecast of Euro (%EUR="Yes")
   'The forecast of Euro as average of Euro area forecasts is stored in "a_feuro"
  ' Hard coded numbers in c(600+!j) for !j=1,..,N1, c(650), c(650) are the place of some coefficient 
  ' estimates stored transitorily in the vector "c" of the workfile.
  
   if %EUR="Yes" then
	smpl @all
	do POOL6.ls
	for !j=10 to N1
	   	if M=1 then
    		scalar sef{!j}= c(600+!j) + c(650)*Fs_{!j}(PFE)
		endif
		if M=2 then
    		scalar sef{!j}= c(600+!j) + c(650)*Fs_{!j}(PFE)+c(660)*xtr_{!j}(PFE)
		endif
		if M=3 then
	    	scalar sef{!j}= c(600+!j) + c(650)*Fs_{!j}(PFE)+c(660)*xmm_{!j}(PFE)
		endif
		if M=4 then
    		scalar sef{!j} = c(600+!j) + c(650)*Fs_{!j}(PFE)+c(660)*xppp_{!j}(PFE)
		endif
	next

	for !j=10 to N1
      		she=she+sef{!j}  
	next
    	a_feuro(!hh)=she/8
     endif

'********************************************************************************************************
'********************************************************************************************************
'
' 6. Forecast evaluation
'
'********************************************************************************************************
' This section constructs prediction errors, Theil's U stats, Clark-West stats for univariate
' and multivariate cases, Max t-stats (maximum of the N1 Clark-West t-statistics) and p-values
' Uncomment line below "RMSPE of RW model " if you want RMSPE of random walk model
'  also stored in the workfile
' Notation:
' pe_ prediction errors  for model (1,...,4) and RW model
' foo_ denote squared prediction errors 
' dhm_ and  dbarm_ are auxiliary variables to construct CW stats below
' svb is an adjusted version of P for CW stats in multivariate case
' gt, gbar, ght, vhs, vhsx etc are auxiliary matrices or series to construct CW stats
' vhats is long run variance estimated using subroutine (see section 1)
' zs is CW stat for each country; zhats=max of z-stats
' p_vals are p-values for max z-stat
' Hard coded number: in "!j>9", actually from the 10th currency on, for countries in Euro area
' NX1=9 (non Euro currencies), N1 *(all), 10 (non Euro currencies+ Euro) for early, full, late 
' sample for CW stats in multivariate case
'
' Setting sample

 if %S="E" then
 smpl %fp+t0    1998q4 
 endif

 if %S="L" then
 smpl 1999q1    @last-k
 endif 
 
 if %S="F" then
 smpl %fp+t0    @last-k
 endif 

'********************************************************************************************************
'  Prediction Errors and U stats (and RMSPE if needed) for j=1,...,N1

 for !j = 1 to N1

   if %S="L" and !j>9 then 'Euro Prediction Error 
   series pe_m1_{!j} = (s_{!j}(+k)-s_{!j})-shateuro_{!hh}
   else
   series pe_m1_{!j} = (s_{!j}(+k)-s_{!j})-shat_{!j}
   endif
      
   series pe_rw{!j} = (s_{!j}(+k)-s_{!j}) 
   series foo_m1=pe_m1_{!j}*pe_m1_{!j} 
   series foo_rw=pe_rw{!j}*pe_rw{!j}  
 
   if %S="L" and !j>9 then
   series dhm1_{!j}=(foo_rw)-(foo_m1)+(shateuro_{!hh}*shateuro_{!hh})
   else
   series dhm1_{!j}=(foo_rw)-(foo_m1)+(shat_{!j}*shat_{!j})
   endif

   series foo_m1_k{!hh}=foo_m1
   series foo_rw_k{!hh}=foo_rw
   series dhm1_{!j}_k{!hh}=dhm1_{!j}   

   scalar dbarm_{!j} = @mean(dhm1_{!j})

   'RMSPE of RW model ('sugg.: uncomment if want RMSPE of random walk model)
   'a_rmspe_rw(!j,!hh) = @mean(foo_rw)^.5   

   'Theil's U stats
   a_ucw_s(!j,!hh)=(@mean(foo_m1)/@mean(foo_rw))^.5
   delete foo_m1 foo_rw 'pe_m1_{!j}

 next

'********************************************************************************************************
' Univariate case: Standard errors and CW stats 

' Setting sample

 if %S="F" then
 smpl %fp+t0    @last-k+1
 scalar P1=P-k+1
 scalar P2=P-2*(k-1)
 endif 

 if %S="E" then
 smpl %fp+t0    1998q4+1 
 scalar P1=P-k+1
 scalar P2=P1
 endif

 if %S="L" then
 smpl 1999q1    @last-k+1
 scalar P1=P-(k-1)
 scalar P2=P-2*(k-1)
 endif 

   for !j=1 to N1
   series ssx_{!j}=0
   	for !g=1 to  k
   	 if %S="L" and !j>9 then 'for Euro std error 
   	  series ssx_{!j}=ssx_{!j}+shateuro_{!hh}(-!g)
	 else
   	  series ssx_{!j}=ssx_{!j}+shat_{!j}(-!g)
	 endif
   	next
   series ssh_{!j}=ssx_{!j}
   series gt_{!j} =2*(s_{!j}-s_{!j}(-1))*ssh_{!j}
   series gt_{!j}_{!hh} =gt_{!j} 
   scalar gbar_{!j}=@mean(gt_{!j})
   series gb_{!j}=gt_{!j}-gbar_{!j}
   scalar vhs_{!j}=(1/P2)*@sumsq(gb_{!j}) 
   vector(N1) vhsx(!j)=vhs_{!j}

   'Univariate Clark-West stats
   scalar tcw_{!j}=@sqrt(P1)*dbarm_{!j}/(@sqrt(vhs_{!j}))
   vector(N1) tcwx(!j)=tcw_{!j}
   delete ssh_{!j} 'ssx_{!j} shat_{!j}	'uncomment/add more auxiliary variables if need delete
   a_ucw_s(!j,!hh+rh)=tcw_{!j}	' storing the CW stats
   next

'********************************************************************************************************
 ' Multivariate case: Standard errors and CW stats 
 ' Estimates long run variance (vhats) for each sample

' Setting sample
 if %S="F" then
 scalar svb=P2
 scalar NX1=9
 endif
 
 if %S="E" then
 scalar svb=P-k+1
 scalar NX1=N1
 endif

 if %S="L" then
 scalar svb=P2
 scalar NX1=10 
 endif

  matrix(NX1,svb) Vbig 
  matrix (NX1,NX1) Vhats

  for !j=1 to NX1 
  vector ght_{!j} =@transpose(@convert(gb_{!j}))
  rowplace(Vbig,ght_{!j},!j)
  delete ght_{!j}  gb_{!j} gt_{!j}
  next

  if %S="F" then
  call gam0(P1,k,NX1,Vbig,Vhats) 
  endif

  if %S="E" then
  call gam0(svb,1,NX1,Vbig,Vhats) 
  endif

  if %S="L" then
  call gam0(P1,k,NX1,Vbig,Vhats)
  endif

  vector vhs = @sqrt(@getmaindiagonal(Vhats))
  sym capDs_1 = @makediagonal(vhs)
  matrix omehats=(@inverse(capDs_1))*Vhats*(@inverse(capDs_1))
  delete capDs_1 Vhats Vbig

'********************************************************************************************************
' Max t-stat (maximum of the N1 Clark-West t-statistics) and p-values

  vector(NX1) zs
  for !j=1 to NX1
   zs(!j)=@sqrt(P1)*dbarm_{!j}/vhs(!j) 	'z-stat(j) for j=1,..., last currency
  'delete dbarm_{!j}
  next
  vector zhats = @cmax(zs)		' max of z-stat(j)
  delete zs  'vhs

   matrix (ND,NX1) zms 
   vector (ND) rzs
   call rdraw(ND,NX1,omehats, zms, rzs)	' generating random draws
   delete omehats zms 

   vector (1) p_vals
   call pvals(ND,zhats, rzs, p_vals)		' computing p-values for max z-stat
   a_ucw_s(N1+2,!hh+rh)=zhats(1)
   a_ucw_s(N1+3,!hh+rh)=p_vals(1)   
   delete zhats rzs

next	' main loop  for each forecast horizon ends here

'********************************************************************************************************
'********************************************************************************************************
'
' 7. Excel Output
'
'********************************************************************************************************
' Saves content specified just before ".write..." in an Excel spreadsheet in the specified directory.
' Change the location and name of directory accordingly.
' In particular, it saves: 
' First, the "a_ucw_s" matrix (U and CW stats) in an Excel sheet with name 
' "stats_[model#]_[#factors].xls" according to values chosen for M (model) and F 
'  (number of factors) specified in section 2 (Main Parameters). That is, it reports only stats for 
' a specific model, a specific number of factors (and a specific sample). Note also that the 
'  "a_ucw_s" matrix is also stored in the workfile.
' The ouput is an arrangemement with 20 rows and 10 columns. 
' The upper left 17x5 submatrix reports the U stats, the upper right 17x5 submatrix reports
'  the CW stats, the lower right 2x5 submatrix reports max t-stat and 
' p-values, the rest of cells are filled with zeros. Note that the 18th row only has zeros. 
' The arrangement is similar to the one used to report results as in emw23.xls.
' Second, it also saves "a_feuro", a 5x1 vector under the name "forecast_[model#]_[#factors].xls"
' that is the vector with the forecast of the Euro if  %EUR="Yes"  in section 2 .

' The following is a dummy loop that helps to index the Excel files per model# and # of factors.

     for !f=F to F
     	for !m=M to M
     		a_ucw_s.write(t=xls) c:\courses\RA\exchanges\emw40b\stats_M{!m}_{!f}F
		if %EUR="Yes" then
		a_feuro.write(t=xls) c:\courses\RA\exchanges\emw40b\forecast_M{!m}_{!f}F
		endif
     	next     
     next

'********************************************************************************************************
'C'EST FINI