
libname home 'wherever the quartile variables sorted data set is located';
%let csvdir= %str(wherever the quartile variables sorted data set is located);
options obs=max;

%let fyr = 1996;

/* import full data and compute all required variables */
%macro add_vars;
  data working_data(drop=ssa_yob gnd_ind change_jobs holding_wages prev_val_fj lagtin);
    set home.quartile_variables_sorted(keep=tin tax_yr payer_tin_w2_max ssa_yob gnd_ind inventor
    wages1996 wages1997 wages1998 wages1999 wages2000 wages2001 wages2002 wages2003 wages2004
    wages2005 wages2006 wages2007 wages2008 wages2009 wages2010 wages2011 wages2012 wages2013 wages2014
    payer_tin1996 payer_tin1997 payer_tin1998 payer_tin1999 payer_tin2000 payer_tin2001 payer_tin2002
    payer_tin2003 payer_tin2004 payer_tin2005 payer_tin2006 payer_tin2007 payer_tin2008
    payer_tin2009 payer_tin2010 payer_tin2011 payer_tin2012 payer_tin2013 payer_tin2014) end=eof;

    /* assign wages and separation in each year */
    wages=.;
    sep=0;
    %do year=&fyr. %to 2014;
      %let nxt = %eval(&year. + 1);
      if tax_yr = &year. then
      do;
        wages = wages&year.;
        if payer_tin_w2_max ne payer_tin&nxt. and tax_yr ne 2014 then sep=1;
        drop wages&year. payer_tin&year.;
      end;
    %end;

    /* DONT NEED: logwages = .;
    if wages ne . and wages > 0 then logwages = log(wages);*/

    /* quadratic age */
    age = tax_yr-ssa_yob;

    /* replace inventor with 0 if missing */
    if inventor=. then inventor=0;

    /* gender */
    male=0;
    if gnd_ind="M" then male=1;
    female=0;
    if gnd_ind="F" then female=1;

    /* previous year's firm */
    lagtin = lag1(payer_tin_w2_max);

    /* previous job's wages */
    change_jobs = (tax_yr>&fyr.)*(payer_tin_w2_max ne lagtin);
    prev_wages = holding_wages;
    holding_wages = lag1(wages);

    if change_jobs=1 and ((lagtin ne . and lagtin ne 0) or prev_val_fj = 1) then prev_wages=holding_wages;
    else if prev_wages=. or  tax_yr=&fyr. then prev_wages=0;
    holding_wages = prev_wages;
    retain holding_wages;

    /* first job */
    if tax_yr = &fyr. then first_job = 1;
    else if change_jobs = 1 and lagtin ne 0 and lagtin ne . then first_job = 0;
    else first_job = prev_val_fj;
    prev_val_fj = first_job;
    retain prev_val_fj;

    /* make variables for regression */
    age_sq = age*age;
    age_cu = age_sq*age;
    age_qu = age_cu*age;
    prev_wages_sq = prev_wages*prev_wages;

    male_age = male*age;
    male_age_sq = male*age_sq;
    male_age_cu = male*age_cu;
    male_age_qu = male*age_qu;

    female_age = female*age;
    female_age_sq = female*age_sq;
    female_age_cu = female*age_cu;
    female_age_qu = female*age_qu;

    inventor_age = inventor*age;
    inventor_age_sq = inventor*age_sq;
    inventor_age_cu = inventor*age_cu;
    inventor_age_qu = inventor*age_qu;

    male_inventor = male*inventor;
    female_inventor = female*inventor;

    /* we do not want values that are zero */
    if payer_tin_w2_max=0 or payer_tin_w2_max=. or wages=0 or wages=. then delete;

  run;
%mend;
%add_vars;

/* import full data and compute all required variables */
%macro predict_vals(var_list, newvar);

  data _null_;
    set ests(keep=Parameter Estimate) end=eof;

      /* declare a hash table */
      if _N_ = 1 then
      do;
        declare Hash estimates ();
          rc = estimates.DefineKey ('tax_yr');
          rc = estimates.DefineData ('tax_yr','Estimate');
          rc = estimates.DefineDone ();
      end;

      /* fill it with values */
      if indexw(Parameter,'tax_yr') then
      do;
        Parameter = translate(Parameter,'','tax_yr');
        Parameter = trim(Parameter);
        tax_yr = 1*Parameter;
        rc = estimates.ref();
      end;

      /*output the hash table */
      if eof then rc = estimates.output(dataset:'estimates');

  run;

  data working_data(drop=rc Estimate);
    set working_data;

    length tax_yr 4 Estimate 8;

    if _N_ = 1 then do;
      declare Hash estimates(dataset:'estimates');
      rc = estimates.DefineKey('tax_yr');
      rc = estimates.DefineData('Estimate');
      rc = estimates.DefineDone();
    end;

    rc = estimates.find();
    &newvar. = &newvar. - Estimate;

  run;

%mend;

/* run the quality index regressions: first baseline then expanded */

/* wrapper for quality regressions */
%macro run_reg(depvar, indepvars, outvar);

  ods trace on;
  proc glm data=working_data;
    class tax_yr;
    model &depvar. = &indepvars. tax_yr / solution;
    output out=working_data predicted=&outvar.;
    ods output ParameterEstimates=ests;
  run;
  ods trace off;
  %predict_vals(var_list=&indepvars., newvar=&outvar.);

%mend;

/* baseline level quality */
%run_reg(depvar=wages, indepvars=age age_sq age_cu age_qu male female inventor male_age male_age_sq male_age_cu male_age_qu female_age female_age_sq female_age_cu female_age_qu inventor_age inventor_age_sq inventor_age_cu inventor_age_qu male_inventor female_inventor, outvar=baseline_quality)

/* expanded level quality */
%run_reg(depvar=wages, indepvars=age age_sq age_cu age_qu male female inventor male_age male_age_sq male_age_cu male_age_qu female_age female_age_sq female_age_cu female_age_qu inventor_age inventor_age_sq inventor_age_cu inventor_age_qu male_inventor female_inventor prev_wages prev_wages_sq first_job, outvar=expanded_quality)

/* baseline sep rates */
%run_reg(depvar=sep, indepvars=age age_sq age_cu age_qu male female inventor male_age male_age_sq male_age_cu male_age_qu female_age female_age_sq female_age_cu female_age_qu inventor_age inventor_age_sq inventor_age_cu inventor_age_qu male_inventor female_inventor, outvar=baseline_sep_prob)

/* expanded sep rates*/
%run_reg(depvar=sep, indepvars=age age_sq age_cu age_qu male female inventor male_age male_age_sq male_age_cu male_age_qu female_age female_age_sq female_age_cu female_age_qu inventor_age inventor_age_sq inventor_age_cu inventor_age_qu male_inventor female_inventor prev_wages prev_wages_sq first_job, outvar=expanded_sep_prob)

/* get the median for each probability into macro variables */
proc sql;

  select med_qual_b, med_qual_e, med_sep_b, med_sep_e
    into       :mqb,       :mqe,      :msb,      :mse
    from (
      select median(baseline_quality)  as  med_qual_b, median(expanded_quality)  as med_qual_e,
             median(baseline_sep_prob) as   med_sep_b, median(expanded_sep_prob) as med_sep_e
        from working_data
    );

  create table qs_bln as
    select payer_tin_w2_max, tax_yr,
           baseline_quality>=&mqb  as abv_med_qual_bln,
           baseline_sep_prob>=&msb as abv_med_sep_bln,
           mean(wages) as w2_bln_qual_sep,
           count(*) as num_bln_qual_sep
      from working_data group by payer_tin_w2_max, tax_yr, abv_med_qual_bln, abv_med_sep_bln;

  create table qs_exp as
    select payer_tin_w2_max, tax_yr,
           expanded_quality>=&mqe  as abv_med_qual_exp,
           expanded_sep_prob>=&mse as abv_med_sep_exp,
           mean(wages) as w2_bln_qual_exp,
           count(*) as num_bln_qual_exp
      from working_data group by payer_tin_w2_max, tax_yr, abv_med_qual_exp, abv_med_sep_exp;

quit;

proc export data=qs_exp outfile="&csvdir./med_expanded_qual_sep_vars.csv" dbms=csv replace;
run;

proc export data=qs_bln outfile="&csvdir./med_baseline_qual_sep_vars.csv" dbms=csv replace;
run;
