/*----------------------------------------------------------------------------
    Date: May 2007
    run:  computing weather normalized volume
    Data: weathernorm2.xls
        : OEB Gas empirical
        : growth rate model
        : growth rates computed using ln(data/datalag)
        : models: gyv = f(ghdd)
        : LHS: gyvrc, gyvoth, gyvcom,gyvres
-----------------------------------------------------------------------------*/

new; et = hsec; 

library dstat, lr, pgraph; 
lrset; graphset;
output file = C:\work\Oebgas\results\wAdj2 reset;
outfile = "C:\\work\\Oebgas\\results\\wAdj2";

screen on;

@Specify temporary files:@

temp_1 = "C:\\work\\Oebgas\\Temp_1.dat"; 
temp_2 = "C:\\work\\Oebgas\\Temp_2.dat"; 
temp_3 = "C:\\work\\Oebgas\\Temp_3.dat"; 
temp_4 = "C:\\work\\Oebgas\\Temp_4.dat"; 

xls_file = "C:\\work\\Oebgas\\weathernorm2.xls";

importf(xls_file,temp_1,0,0);  

? "  ";	
	? "********************************************************************************"; 
	? "  ";	
	? " Date: ";; ? Datestr(0);;
    
    ? " ****    ";; ? "REGRESSION FOR WEATHER ADJUSTING VOLUME DATA";;
    	? "   ****   ";; ? "Time: ";; ? Timestr(0);
	? "  ";
	? "         		OUTPUT FILE:";; ? outfile;
	? "  ";	
    ? "         		    DATA FILE:";; ? xls_file;	
	? " ";     
    ? "********************************************************************************"; 
 
dataloop ^Temp_1 ^Temp_2; 

Select Year >= 1994 and Year <= 2005; 

delete bad == 1;

keep id year bad hdd fhdd yvrc yvoth yvres yvcom;

endata;

open f1 = ^Temp_2;
data = readr(f1,rowsf(f1));

@ Set up here to be general enough for unbalanced panel: @

id_now = Data[.,1];
grpnums = id_now[1,1]; j = 2;
do until j == rows(id_now); 
  if id_now[j,1] /= id_now[j-1,1]; grpnums = grpnums|id_now[j,1]; endif;
j = j + 1; 
endo;
dum = id_now[1:rows(id_now),.] .== grpnums'; T_period = sumc(dum); 
N_firms = rows(T_period); 

@ lag variables and divide original vars by lagged values @
@ Log this ratio to get growth rates of vars    @

{ vnames,indx } = indices(Temp_2,0);
let var_name = hdd;
start = indcv(var_name,vnames);

colsff1 = colsf(f1);
datal= lag1(data ); 

data2 = data[.,1:start-1]~ln(data[.,start:colsff1]./datal[.,start:colsff1]);

f1=close(f1);
open f1 = ^Temp_2 for update; 
writer(f1,data);
f1 = close(f1);

@ 1. create a vector with sequenced values for each cross section
  2. append to data2
  3. use to delete first row of each cross section - i.e. drop first obs of
     each firm that has been lagged @
  
id1 = zeros(rows(data2),1); 
lb = 1; ub = 0; firm = 1;
do until firm > N_firms;
ub = ub + T_period[firm,1];

id1[lb:ub,1] = seqa(1,1,rows(id1[lb:ub,1]));

firm = firm+1; 
lb = ub+1; 
endo;

vnames = getname(Temp_2);
vnames2 = vnames|"rind";                        @ rind is for row index @
create h6 = ^temp_3 with ^vnames2,0,8; 
if h6 == -1; 
errorlog "Can't create data set Temp_3"; end; 
endif; 

data2 = data2~id1;

y = writer(h6,data2); 
h6 = close(h6); 

data2 = delif(data2,(data2[.,indcv("rind",vnames2)] .== 1)); 

@***********************************************************************************@

rhs = { fhdd };

__con = 0;
_olsres = 1;
__altnam = { gfhdd,gyvrc };
x = ones(rows(data2),1)~data2[.,indcv(rhs,vnames2)];
y = data2[.,indcv("yvrc",vnames2)];

? "                           OLS REGRESSION "; ? "  ";
{ vnam, m, byvrc, stdbols, vcols, seols, sigols, cx, rsq, eols, dwstat } = ols(0,y,x);
? "  "; ? "  ";

@***********************************************************************************@

rhs = { fhdd }; 

__con = 0;
_olsres = 1;
__altnam = { gfhdd,gyvoth }; 
x = ones(rows(data2),1)~data2[.,indcv(rhs,vnames2)];
y = data2[.,indcv("yvoth",vnames2)];

? "                           OLS REGRESSION "; ? "  ";
{ vnam, m, byvoth, stdbols, vcols, seols, sigols, cx, rsq, eols, dwstat } = ols(0,y,x);
? "  "; ? "  ";

@***********************************************************************************@

rhs = { fhdd }; 

__con = 0;
_olsres = 1;
__altnam = { gfhdd,gyvres }; 
x = ones(rows(data2),1)~data2[.,indcv(rhs,vnames2)];
y = data2[.,indcv("yvres",vnames2)];

? "                           OLS REGRESSION "; ? "  ";
{ vnam, m, byvres, stdbols, vcols, seols, sigols, cx, rsq, eols, dwstat } = ols(0,y,x);
? "  "; ? "  ";

@***********************************************************************************@

rhs = { fhdd }; 

__con = 0;
_olsres = 1;
__altnam = { gfhdd,gyvcom }; 
x = ones(rows(data2),1)~data2[.,indcv(rhs,vnames2)];
y = data2[.,indcv("yvcom",vnames2)];

? "                           OLS REGRESSION "; ? "  ";
{ vnam, m, byvcom, stdbols, vcols, seols, sigols, cx, rsq, eols, dwstat } = ols(0,y,x);
? "  "; ? "  ";

@***********************************************************************************@


weatherA = zeros(rows(data),6); 

lb = 1; ub = 0; firm = 1;
do until firm > N_firms;
ub = ub + T_period[firm,1];

yvrc=data[.,indcv("yvrc",vnames2)];
yvoth=data[.,indcv("yvoth",vnames2)];
yvres=data[.,indcv("yvres",vnames2)];
yvcom=data[.,indcv("yvcom",vnames2)];
fhdd=data[.,indcv("fhdd",vnames2)];

weatherA[lb:ub,1] = data[lb:ub,indcv("id",vnames2)];
weatherA[lb:ub,2] = data[lb:ub,indcv("year",vnames2)];
weatherA[lb:ub,3] = exp(ln(yvrc[lb:ub,.])+byvrc[2,.]*(ln(meanc(fhdd[lb:ub,.])./fhdd[lb:ub,.])));
weatherA[lb:ub,4] = exp(ln(yvoth[lb:ub,.])+byvoth[2,.]*(ln(meanc(fhdd[lb:ub,.])./fhdd[lb:ub,.])));
weatherA[lb:ub,5] = exp(ln(yvres[lb:ub,.])+byvres[2,.]*(ln(meanc(fhdd[lb:ub,.])./fhdd[lb:ub,.])));
weatherA[lb:ub,6] = exp(ln(yvcom[lb:ub,.])+byvcom[2,.]*(ln(meanc(fhdd[lb:ub,.])./fhdd[lb:ub,.])));

firm = firm+1; 
lb = ub+1; 
endo;

wAdj2 = "C:\\work\\OEBGAS\\wAdj2.xls";
name ="id"$~"year"$~"Ayvrc"$~"Ayvoth"$~"Ayvres"$~"Ayvcom";
call xlswritesa(name,wAdj2,"a1",1,"");
call xlswritem(weatherA,wAdj2,"a2",1,0);

closeall; Output off;

/*----------------------------------- END of PROGRAM ------------------------------*/