/*
              FILE: persmooK_RS.pgm
              This version:  September 21, 2000.
*/
new ; cls ;
clear test0,test1,test2,test3,test4 ;
ridge = .10 ; @weight of ridge in GLS @

off=1; pc=0 ; bos = 0 ;

if pc==1 ; load path=c:\aba\fdic\data\fmt ; endif ;
if bos==1 ; load path=m:\fdic\data\fmt ; endif ;
if off==1 ; load path=d:\fdic2001\data\fmt ; endif ;


/* CONTROL VARIABLES */
        _nondur = 1 ;                 /* nondurable consumption */
              _logs = 1 ;                 /* ln variables */
            
ngroups=1 ;

if pc==1 ;output file=c:\aba\fdic\sept00\smoo_rs.out         reset; endif;
if bos==1;output file=m:\fdic\sept00\smoo_rs.out             reset; endif;
if off==1;output file=d:\fdic2001\sept00\smoo_rs.out reset; endif;


"program file: persmooK_RS.pgm " ;
"date and time : " ;;  datestr(date) ;;" ";;timestr(time) ;
"relative boosting of diagonal is: " ridge ;

	"LOGS ON " ;
 " " ;

/* LOAD DATA */
/* All fmt data files in milions */

stateno = seqa(1,1,51) ;
stateno = stateno[1:8 10:51,.] ;

load pop6396 ;
pops=pop6396[1:8 10:51,.] ;
pops = pops[.,1:33] ;
clear pop6396 ;

load cpi6396 ;
cpi_vec = cpi6396[1:33] ;
cpi=ones(rows(pops),1).*.cpi_vec' ;  /* 63 - 95 */
clear cpi6396 ;

load dpi6396 ;
dpinc = dpi6396[1:8 10:51,.] ;
dpinc = dpinc[.,1:33] ;
clear dpi6396 ;

dpi_agg = sumc(dpinc) ;
dpinc = dpinc./cpi ;                  
dpi_agg = dpi_agg./cpi_vec ;


"NONDURABLE RETAIL SALES" ;
load ndur6095 ;
ndur6395 = ndur6095[.,4:36] ;
ndur = ndur6395[1:8 10:51,.]  ;
clear ndur6395 ;
sale = ndur./cpi ;
perc = sumc(ndur) ;
perc = perc./cpi_vec ;                  
"AGG CONSUMPTION IS SUM OF NONDURABLE RETAIL SALES" ;


/* PER CAPITA  */
dpi = dpinc./pops ;
sale = sale./pops ;
perc = perc./sumc(pops) ;
dpi_agg = dpi_agg./sumc(pops) ;
dpi_agg = ones(rows(pops),1).*.dpi_agg' ;  
perc = ones(rows(pops),1).*.perc' ;  

dpi_l = ln(dpi) ;
sale_l = ln(sale) ;
perc_l = ln(perc) ;
dpi_aggl = ln(dpi_agg) ;

/* DIFFERENCES */
p=1 ; "length of difference : ";; p ;

capT =  cols(dpi_l) ; 
capN = rows(dpi_l) ;

dpi_d    = delta_bent(dpi_l,p,capT) ;
sale_d  = delta_bent(sale_l,p,capT) ;
perc_d  = delta_bent(perc_l,p,capT) ;
dpi_aggd  = delta_bent(dpi_aggl,p,capT) ;

/* LAGS */
kmax = 1 ;   	   /* max number of lags */

dpi_t  = dpi_d' ;
dpi_1 = lagn(dpi_d',1) ;
dpi_2 = lagn(dpi_d',2) ;
dpi_t  = trimr(dpi_t, kmax , 0)' ;
dpi_1 = trimr(dpi_1, kmax ,0)' ;
dpi_2 = trimr(dpi_2, kmax ,0)' ;

sale_t  = sale_d' ;
sale_1 = lagn(sale_d',1) ;
sale_2 = lagn(sale_d', 2) ;
sale_t = trimr(sale_t, kmax ,0)' ;
sale_1 = trimr(sale_1, kmax ,0)' ;
sale_2 = trimr(sale_2, kmax ,0)' ;

dpi_aggt  = dpi_aggd' ;
dpi_agg1 = lagn(dpi_aggd',1) ;
dpi_agg2 = lagn(dpi_aggd',2) ;
dpi_aggt = trimr(dpi_aggt, kmax ,0)' ;
dpi_agg1 = trimr(dpi_agg1, kmax ,0)' ;
dpi_agg2 = trimr(dpi_agg2, kmax ,0)' ;

perc_t  = perc_d' ;
perc_1 = lagn(perc_d',1) ;
perc_2 = lagn(perc_d',2) ;
perc_t = trimr(perc_t, kmax ,0)' ;
perc_1 = trimr(perc_1, kmax ,0)' ;
perc_2 = trimr(perc_2, kmax ,0)' ;

/** OUTPUT /
format /m1 /rd 5,2 ;

N = rows(dpi_t) ; T = cols(dpi_t) ; one_w=ones(N,T) ;
"fixfax with nationwide average" ;
"sample: 1963-95";
" " ;

_timfix = 1 ; 
_crsfix = 0 ; 
if (_timfix==1); "time fixed eff"; endif ;
if (_crsfix==1); "cross section fixed eff"; endif ;

{dpi_t_ft,dpi_1_ft} = fixfax(dpi_t,dpi_1,one_w) ;
{dpi_atft,dpi_a1ft}= fixfax(dpi_t-dpi_aggt,dpi_1-dpi_agg1,one_w) ;
{saleatft,perc_tft} = fixfax(sale_t-perc_t,perc_t,one_w) ;
{sale_tft,perc_tft} = fixfax(sale_t,perc_t,one_w) ;

_timfix = 0 ; 
_crsfix = 1 ; 
if (_timfix==1); "time fixed eff"; endif ;
if (_crsfix==1); "cross section fixed eff"; endif ;

{dpi_t_fx,dpi_1_fx} = fixfax(dpi_t,dpi_1,one_w) ;
{dpi_atfx,dpi_a1fx}= fixfax(dpi_t-dpi_aggt,dpi_1-dpi_agg1,one_w) ;
{sale_tfx,perc_tfx} = fixfax(sale_t,perc_t,one_w) ;
{saleatfx,perc_tfx} = fixfax(sale_t-perc_t,perc_t,one_w) ;

_timfix = 1 ; 
_crsfix = 0 ; 

{dpi_t_fxt,dpi_1_fxt} = fixfax(dpi_t_fx,dpi_1_fx,one_w) ;
{sale_tfxt,perc_tfxt} = fixfax(sale_tfx,perc_tfx,one_w) ;
{saleatfxt,perc_tfxt} = fixfax(saleatfx,perc_tfx,one_w) ;

format /m1 /rd 5,2 ;
"  BENCHMARK **" ;
" " ;
"Regression of C_it on constant, Y_it w NO fix eff" ;
" " ;
{gls,glsstdv} = pest(one_w~dpi_t,sale_t) ;     /* pest(x,y) */
" " ;
"Coefficients:        " ;
"b_gls  :     " ;;  gls' ;
"SE_gls :     " ;; glsstdv' ;
"t-test :     " ;;  (gls./glsstdv)' ;
" " ;

" * RISK SHARING ***" ;
" " ;
"Regression of C_it on Y_it with TIME fixed effects and no cross f.ex." ;
" " ;
{gls,glsstdv} = pest(dpi_t_ft,sale_tft) ;     /* pest(x,y) */
" " ;
"Coefficients:        " ;
"b_gls  :     " ;;  gls' ;
"SE_gls :     " ;; glsstdv' ;
"t-test :     " ;;  (gls./glsstdv)' ;
" " ;
"Regression of (C_it - agg C_t) on constant, (Y_it - agg Y_t) w NO fix eff" ;
" " ;
{gls,glsstdv} = pest(one_w~(dpi_t-dpi_aggt),(sale_t-perc_t)) ;     /* pest(x,y) */
" " ;
"Coefficients:        " ;
"b_gls  :     " ;;  gls' ;
"SE_gls :     " ;; glsstdv' ;
"t-test :     " ;;  (gls./glsstdv)' ;
" " ;
" " ;
"Regression on LAGGED income With time fixed effects " ;
{gls,glsstdv} = pest(dpi_1_ft,sale_tft) ;           /* olswin(x,y) */
" " ;
"Coefficients:        " ;
"b_gls  :     " ;;  gls' ;
"SE_gls :     " ;; glsstdv' ;
"t-test :     " ;;  (gls./glsstdv)' ;
" " ;
" " ;

"With  Aggregate income AND agg consumption subtracted  " ;
"delta (c_it - C_it) = a_0 + b1 delta (y_it-1 - Y_t-1) + error" ;
" " ;

{gls,glsstdv} = pest(one_w~(dpi_1-dpi_agg1),(sale_t-perc_t)) ;     /* pest(x,y) */
" " ;
"Coefficients:        " ;
"b_gls  :     " ;;  gls' ;
"SE_gls :     " ;; glsstdv' ;
"t-test :     " ;;  (gls./glsstdv)' ;




" * SMOOTHNESS REGRESSIONS ***" ;
" " ;
"Regression of (C_it - agg C_t) on (Y_it - agg Y_t) w cross fix eff" ;
" " ;
{gls,glsstdv} = pest(dpi_atfx,saleatfx) ;     /* pest(x,y) */
" " ;
"Coefficients:        " ;
"b_gls  :     " ;;  gls' ;
"SE_gls :     " ;; glsstdv' ;
"t-test :     " ;;  (gls./glsstdv)' ;
" " ;


"Regression of C_it on Y_it and both cross and  time fixed effect " ;
{gls,glsstdv} = pest(dpi_t_fxt,sale_tfxt) ;     /* pest(x,y) */
" " ;
"Coefficients:        " ;
"b_gls  :     " ;;  gls' ;
"SE_gls :     " ;; glsstdv' ;
"t-test :     " ;;  (gls./glsstdv)' ;
" " ;

output off ;

/***
    PROCEDURES
***/


/**** PROCEDURES ****/
proc(1)=delta_bent(xxx,diflen,lastyear) ;
local s,dx,i,k ;
k = floor( (lastyear-1)/diflen )  ;
s=0 ; i=0 ;
dx=zeros(rows(xxx),k) ;
do while s<(lastyear-diflen ) ;
   dx[.,k-i]= xxx[.,lastyear-s] - xxx[.,lastyear-s-diflen] ;
   s=s+diflen ; i=i+1 ;
endo ;
retp(dx) ;
endp ;


/** fixed effects removal***/
proc(2)=fixfax(y,x1,het_w) ;
local hetvec,yvec,xvec,TT,NN,nobs,x1vec,yNs,x1Ns,yTs,x1Ts ;
local olsest,xx,madest,stdev,resid,sighat,olsdv ;

/ Note: No diff in treatment of x and y /
NN=rows(y)  ;
TT=cols(y) ;
Nobs=TT*NN ;

/ heteroskedasticity *
Heteroskedasticity weigths (inverse of the std devs)   /

if _timfix==1 ;
 yNs = sumc(y.*het_w)./sumc(het_w^2) ;  @yNs: TTx1 - for each t sum over n@
 y = y - het_w.*(ones(NN,1).*.yNs') ;     @time fixed effect@

 x1Ns = sumc(x1.*het_w)./sumc(het_w^2) ;  @yNs: TTx1 - for each t sum over n@
 x1 = x1 - het_w.*(ones(NN,1).*.x1Ns') ;     @time fixed effect@
endif ;

if _crsfix==1 ;
 yTs = sumc((y.*het_w)')./sumc(het_w'^2);   @yTs: NNx1 - for each n sum over t@
 y = y - het_w.*(yTs.*.ones(1,TT)) ;     @cross fixed effect@

 x1Ts = sumc((x1.*het_w)')./sumc(het_w'^2); @yTs: NNx1 - for each n sum over t@
 x1 = x1 - het_w.*(x1Ts.*.ones(1,TT)) ;     @cross fixed effect@
endif ;

retp(y,x1) ;
endp ;



proc(2) = gamloop(resNT) ; 
local NN,ii,jj,gamt,res,varians ;

NN = rows(resNT) ;

/* construct Gammat for spatial  correlation  */
/* note that gammat is normalized to be correlations */
/* so the data needs to be transformed correspondingly */
Gamt = eye(NN) ;
varians = ones(NN,1) ;
ii=1 ;
do while  ii<(NN+1)  ;
      varians[ii] = meanc(resNT[ii,.]'.^2) ;
         ii=ii+1 ;
endo ;

ii=1 ;
do while  ii<(NN+1)  ;
  jj=ii+1 ;
  do while ( jj < (NN+1) );
     Gamt[ii,jj] = ( meanc(resNT[ii,.]'.*resNT[jj,.]')) ;
 	 gamt[ii,jj]=gamt[ii,jj]/( sqrt(varians[ii])*sqrt(varians[jj]) ) ;
     Gamt[jj,ii] = Gamt[ii,jj] ;
     jj = jj+1 ;
  endo ;
     Gamt[ii,ii] = 1 ;
  ii = ii + 1 ;
endo ;

retp(gamt,varians) ;
endp ;




/** panel estimation ***/
proc(2)=pest(xxreg,yy) ;
local one_w ,TT,NN,nobs,kx,i,j,konstant,t_index,temp,idmat,olscoef ;
local y1,x1,y2,x2,x3,xvec,yvec,ee1,eei,xall,olsstd,ee2,tstat ;
local x1vec,y1vec,y2vec,x2vec,res2,sighat2,x1old,x2old;
local res1,resNT1,vari1,varihat1,gls1,glsdv1,glsi,glsdvi,glsdv ;
local rho,xvxmat,xvymat,glsvar,glsvar_inv,y2vec_t,x2vec_t,temp2 ;
local Gammat,Gam_square,Pmat, x3vec,y3vec,x1new,x2new,ee3 ;
local y_vec,x1_vec,x2_vec,ynew,xnew,yold,xold,resTN1,het_w,time_w ;
local xomegax,xomegay,varians,xbig,indexbig,xnn,ynn,xx,xallvec,chisq ;
NN=rows(yy)  ;
TT=cols(yy) ;
Nobs=TT*NN ;

kx = cols(xxreg)./TT ;
one_w = ones(NN,TT) ;
		
/* **** regressors */

/* OLS regressions  */
      y1vec = reshape(yy,nobs,1) ;
i=1 ;
  xallvec = reshape(xxreg[.,1:i*TT],nobs,1) ;
do until i >= kx ;
    i=i+1 ;
	xallvec = xallvec~reshape(xxreg[.,(TT*(i-1))+1:i*TT],nobs,1)  ;
endo ;

olscoef= inv(xallvec'xallvec)*xallvec'y1vec  ;
ee1=y1vec - xallvec*olscoef ;
resNT1=reshape(ee1,NN,TT) ;

"olscoef " olscoef' ; 


/* Covariance matrix for entire (N*T) system is    */
/* Omega = I_T.*.Gamma    : block-diagonal           */                                   
/* Hence cross correlation matrix is NxN                */            

{Gammat,varians} = gamloop(resNT1) ;

/* rescale x and y so resid has variance 1 */

    vari1 = 1/sqrt( varians ) ;
    varihat1 = vari1.*.(ones(1,TT)) ;
    het_w = reshape(varihat1,nobs,1) ;

      y1vec = y1vec.*het_w ;

i=1 ;
      x2 = xallvec[.,1].*het_w ;
do until i >= kx ;
    i=i+1 ;
	temp = xallvec[.,i].*het_w ;
	x2 = x2~temp  ;
endo ;
xallvec=x2 ;

gammat= ridge*eye(NN)+ (1-ridge)*gammat ; 
Gammat=inv(gammat) ;
@pmat = chol(Gammat) ;@

/* rescale x and y so resid has variance 1 */


y1=reshape(y1vec,NN,TT) ;

x1=reshape(xallvec[.,1],NN,TT) ; 
xall=x1 ;
i=1 ;
do until i >= kx ;
    i=i+1 ;
	temp = reshape(xallvec[.,i],NN,TT) ;
	xall = xall~temp  ;
endo ;

/* X Omega-invers X and X Omega-invers Y  group-by-group */

t_index = 1 ;
xnn=xall[.,t_index] ;
i=1 ;
do until i >= kx ;
    i=i+1 ;
	xnn = xnn~xall[.,(i-1)*TT+t_index]  ;
endo ;
ynn= y1[.,t_index] ;
xomegax = xnn'gammat*xnn ;
xomegay = xnn'gammat*ynn ;
do until t_index > (TT-1) ;
    t_index = t_index+1 ;
    i=1 ;
    xnn=xall[.,t_index] ;
    do until i >= kx ;
        i=i+1 ;
	    xnn = xnn~xall[.,(i-1)*TT+t_index]  ;
    endo ;
    ynn=y1[.,t_index] ;
    xomegax = xomegax+xnn'gammat*xnn ;
    xomegay = xomegay+xnn'gammat*ynn ;
endo ;

gls1=inv(xomegax)*xomegay ;


glsdv1=sqrt(diag(inv(xomegax))) ;


retp(gls1,glsdv1) ;
endp ;