options linesize=80;
options pagesize=60;

data coalgas;
  infile "./coalgaspq1993.prn";
  input state $ coalp coalq gasp gasq
        pop income cdd hdd ctop5 cadj
        gtop5 gadj;

  if state eq "DC" then delete;
  if state eq "AK" then delete;
  if state eq "HI" then delete;


  /* logging price, quantity, etc */
  ln_coalp = log(coalp);
  ln_coalq = log(coalq);
  ln_gasp  = log(gasp);
  ln_gasq  = log(gasq);
  ln_cdd   = log(cdd);
  ln_hdd   = log(hdd);
  ln_Y     = log(income);
  ln_pop   = log(pop);

  drop coalp coalq gasp gasq cdd hdd income pop;

proc print noobs;
  var state gtop5 gadj;

proc means;
proc corr cov;

/* Let's try to estimate the demand curve by OLS */
proc reg;
  model ln_gasp = ln_gasq;
  model ln_gasp = ln_gasq ln_pop ln_Y ln_cdd ln_hdd;
title 'Demand curve by OLS';

/* Let's try to estimate the supply curve by OLS */
proc reg;
  model ln_gasp = ln_gasq;
  model ln_gasp = ln_gasq gtop5 gadj;
title "Supply curve by OLS";

/* Let's try the reduced form */
proc reg;
  model ln_gasp = gtop5 gadj ln_Y ln_pop ln_cdd ln_hdd;
  model ln_gasq = gtop5 gadj ln_Y ln_pop ln_cdd ln_hdd;
title "Reduced form by OLS";


/* Let's try D-curve with both instruments */
proc syslin 2SLS;
  endogenous ln_gasq;
  instruments ln_Y ln_pop ln_cdd gadj gtop5;
  demand: model ln_gasp = ln_gasq ln_Y ln_pop ln_cdd ln_hdd;
  title "Demand curve by IV";

/* Let's try S-curve with all instruments */
proc syslin 2SLS;
  endogenous ln_gasq;
  instruments gadj gtop5 ln_cdd ln_hdd ln_Y ln_pop;
  supply: model ln_gasp = ln_gasq gtop5 gadj;
  title "Supply curve by IV";

/* Let's try 3SLS */
proc syslin 3SLS;
  endogenous ln_gasq ln_gasp;
  instruments gadj gtop5 ln_cdd ln_hdd ln_Y ln_pop;
  supply: model ln_gasp = ln_gasq gtop5 gadj;
  demand: model ln_gasp = ln_gasq ln_Y ln_pop ln_cdd ln_hdd;
  title "Demand and Supply curves by IV";


/* Now, let's include coal prices in gas demand with instruments */

proc syslin 2SLS;
  endogenous ln_gasq ln_gasp ln_coalp;
  instruments gadj gtop5 cadj ctop5 ln_cdd ln_hdd ln_Y ln_pop;
  demand: model ln_gasp = ln_gasq ln_coalp ln_Y ln_pop ln_cdd ln_hdd;
  title "Demand curve including coal price by IV";

proc syslin 3SLS;
  endogenous ln_gasq ln_gasp ln_coalp;
  instruments gadj gtop5 cadj ctop5 ln_cdd ln_hdd ln_Y ln_pop;
  supply: model ln_gasp = ln_gasq gtop5 gadj;
  demand: model ln_gasp = ln_gasq ln_coalp ln_Y ln_pop ln_cdd ln_hdd;
  title "Demand and suppluy curves including coal price by IV";

run;