'==================================================================================================
'==================================================================================================
'==================================================================================================

subroutine local fernandez (series y_lf_ser, series y_hf_ser, series x_hf_ser, scalar ta, scalar s)

' ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
'	PURPOSE:	Temporal disaggregation using Fernandez method
' ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
'	SYNTAX:		series y_hf_ser
'					call fernandez (y_lf_ser, y_hf_ser, x_hf_ser, ta, s)
' ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
'	OUTPUT:		y_hf_ser	high frequency estimate series name
' ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
'	INPUT:		y_lf_ser	low frequency series name, series must be disaggregated
'								equaly on subperiods
'					x_hf_ser	high frequency related series name
'					ta:			type of temporal disaggregation
'					ta=1		sum (flow)
'					ta=2		average (index)
'					ta=3 		last  element (stock) - interpolation
'					ta=4 		first  element (stock) - interpolation
'					s:			number of high frequency points 
'								for each low frequency data points (freq. conversion)
'					s=4		annual to quarterly
'					s=12		annual to monthly
'					s=3		quarterly to monthly
' ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
'	LIBRARY:		d:\Eviews_library\agg_help.prg: AGG()
' ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
'	writen by:
'	Slawomir Dudek
'	MOF 
' ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
'	REFERENCE:	1) Quilis, E.M (2002) "A MATLAB library of temporal disaggregation methods: summary"
'						    Instituto Nacional de Estadistica, Spain
'						2) Fernandez, R.B. (1981) "Methodological note on the estimation of time series"
'						    Review of Economic and Statistic, vol. 63, n. 3, p. 471-478.
' ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

'	Asiggning series to vectors, only not NA observations
'	Determination of forward extrapolation period: !n_up
'	Determination of backward extrapolation period: !n_dn

smpl @all
series t=@trend+1

sample NA_smpl_y @all if y_lf_ser <> NA
sample NA_smpl_x @all if x_hf_ser <> NA

smpl NA_smpl_x
!obs_f_x=@min(t)
!obs_l_x=@max(t)

smpl NA_smpl_y
!obs_f_y=@min(t)
!obs_l_y=@max(t)

!n_up=!obs_l_x-!obs_l_y
!n_dn=!obs_f_y-!obs_f_x

' ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

'	Control checking wether related series cover aggregate sample

if (!n_up<0 OR !n_dn<0) then

statusline  ERROR!!!  Fernandez subroutine() !!!*** Related series do not cover aggregate sample ***!!!
stop

else

' ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

'	Assigning series to matrix objects, only for NA observations

stomna(y_lf_ser, y_lf_temp, NA_smpl_y)
stomna(x_hf_ser, x_hf_temp, NA_smpl_x)

' ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

'	Size of the problem

!n_hf=@rows(y_lf_temp)		
!n_lf=!n_hf/s									' Size of low frequency input


!n_hf_x=@rows(x_hf_temp)					' Size of high frequency input (number of observations)
!p_hf_x_temp=@columns(x_hf_temp)		' Size of high frequency input (number of variables without intercept)
!p_hf_x=!p_hf_x_temp+1					' Size of high frequency input (number of variables with intercept)

' ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

'	Aggregation of low frequency series equaly distributed in subperiods

matrix AGG1
call aggreg (AGG1, 1, !n_lf, s, 0, 0)
matrix y_lf=AGG1*y_lf_temp
delete AGG1

' ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

'	Generation of aggregation matrix AGG [!n_lf x !n_lf*s] for aggregation in procedure

matrix AGG
call aggreg (AGG, ta, !n_lf, s, !n_up, !n_dn)

' ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

'	Preparing the X matrix: including intercept

matrix(!n_hf_x, !p_hf_x) x_hf=1
matplace(x_hf, x_hf_temp , 1, 1)

' ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

'	Temporal aggreagation of the indicator

matrix x_lf=AGG*x_hf

' ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

'	Estimation 

matrix I=@identity(!n_hf_x)
matrix w=I

matrix(!n_hf_x, !n_hf_x) LL=0
for !i=2 to !n_hf_x
	LL(!i, !i-1)=-1		'Auxiliary matrix useful to simplify computations
next

matrix AUX=(I+LL)
matrix w=@inverse(@transpose(AUX)*AUX)										' High frequency VCV matrix (without sigma_a)
matrix WW=AGG*w*@transpose(AGG)											' Low frequency VCV matrix (without sigma_a)
matrix Wi=@inverse(WW)
matrix beta=@inverse(@transpose(x_lf)*Wi*x_lf)*@transpose(x_lf)*Wi*y_lf		' beta ML estimator
matrix u_lf=y_lf-x_lf*beta															' Low frequency residuals
matrix scp=@transpose(u_lf)*Wi*u_lf												' Weighted least squere
matrix sigma_a=scp*(1/(!n_lf-!p_hf_x))												' sigma_a ML estimator
matrix FM=w*@transpose(AGG)*Wi												' Filtering matrix
matrix u_hf=FM*u_lf																	' High frequency residuals

' ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

'	Temporaly disaggregated time series

matrix y_hf=x_hf*beta+u_hf

' ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

'	Information criteria
'	Note: (!p_hf_x+1) is expanded to include the innovational parameter

scalar aic=@log(sigma_a(1,1))+2*(!p_hf_x+1)/!n_lf
scalar bic=@log(sigma_a(1,1))+@log(!n_lf)*(!p_hf_x+1)/!n_lf

' ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

'	Computation of VCV matrix of high frequency estimates

matrix sigma_beta=sigma_a(1,1)*@inverse(@transpose(x_lf)*Wi*x_lf)

matrix VCV_y=sigma_a(1,1)*(@identity(!n_hf_x)-FM*AGG)*w+(x_hf-FM*x_lf)*sigma_beta*@transpose(x_hf-FM*x_lf)

matrix d_y = @sqrt(@getmaindiagonal(VCV_y))		' Std. dev. of high frequency estimate
matrix y_li=y_hf-d_y										' Lower lim. of high frequency estimate
matrix y_ls=y_hf+d_y									' Upper lim. of high frequency estimate

' ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

'	Assigning vector to series

mtos(y_hf, y_hf_ser, NA_smpl_x)

endif

smpl @all

endsub

'==================================================================================================
'==================================================================================================
'==================================================================================================

