data sasuser.val;
input value r1 r2;
valuen = (r1+r2)/2.0;
keep value valuen;
cards;
01 0 10000
02 10000 14999
03 15000 19999
04 20000 24999
05 25000 29999
06 30000 34999
07 35000 39999
08 40000 49999
09 50000 59999
10 60000 69999
11 70000 79999
12 80000 89999
13 90000 99999
14 100000 124999
15 125000 149999
16 150000 174999
17 175000 199999
18 200000 249999
19 250000 299999
20 300000 399999
21 400000 499999
22 500000 749999
23 750000 999999
24 1000000 1000000
;
run;
proc sort data=sasuser.val;
by value;
run;

data state; set sasuser.NJ;
keep value puma1 hht finc hinc hweight persons yrbuilt grent sex age
                   rooms smocapi grapi bldgsz VEHICL smi lhinc;
if unittype=0;
if smocapi > 0 then smi = 1;
else if grapi > 0 then smi = 0;
else delete;
lhinc = log(1+hinc);

run;
proc sort data=state;
by value;

data state; merge state sasuser.val;
by value;
run;

proc sort data=state; by puma1;
data faminc; set state;
if hht=1 or hht=2 or hht=3;
run;

proc means median data=faminc noprint;
 var finc;
 weight hweight;
 by puma1;
 output out=c median=med;
run;

data c;
  set c;
  keep puma1 med;
  run;

data state;
  merge state c;
  by puma1;

data state; set state;
m4 = med*0.8;
/* if m4 > 50200 then m4 = 50200; */
hlinc = 0;
 units=1;
 if persons <= 4 and hinc < m4*(1-(4-persons)*0.1) then hlinc = 1;
 if persons > 4 and hinc < m4*(1+(persons-4)*0.08) then hlinc = 1;
/***************
* This calculates the percentage of households
* below the 80pct of median
******************/
/* Crowded */
uncrowd = 0;
if persons <= rooms then uncrowd = 1;
/* new construction */
new = 0;
if yrbuilt=1 or yrbuilt=2 or yrbuilt=3 then new =1;
/* Affordable */
Aff = 0;
if smocapi > 0 then grapi=.;
if grapi > 0 then smocapi=.;
if smocapi ne . and  smocapi <= 30 then Aff=1;
if grapi ne . and grapi <= 30 then Aff=1;
xnum = Aff*uncrowd*new*hlinc;
enum = Aff*new*hlinc;
cnum = uncrowd*new*hlinc;
denom = new;
denom2 = new*hlinc;
tot = 1;
run;

proc means data=state noprint;
var denom denom2 xnum enum cnum hlinc tot;
weight hweight;
output out = res sum= sdenom sdenom2 snum senum scnum slinc stot;
run;

data res; set res;
rat = snum/sdenom;
rat2 = snum/sdenom2;
rat3 = slinc/stot;
state = 'NJ';
keep state snum senum scnum sdenom sdenom2 slinc stot rat rat2 rat3;
run;

proc print;
run;

data sasuser.owner;  set state; if smocapi > 0 and smocapi <= 99 and hlinc=1;
lhinc = log(5000+hinc); 
sm = log( smocapi/(100-smocapi)); 
run;

data sasuser.rental; set state; if grapi>0 and grapi <= 99 and hlinc=1;
lhinc = log(5000+hinc);
gra = log( (grapi+5)/(110-grapi)); 
run;

 proc reg data= sasuser.rental;
 weight hweight;
model gra= persons new uncrowd sex age
               rooms  bldgsz VEHICL lhinc;
output out=b cookd=cd p=p r =r;
run;

proc gplot data =b;
plot cd*r r*p;
run;

data sasuser.nym;
infile 'c:\documents and settings\a\My Documents\courses\486n586\nym.txt';
input nn ST CO EXPEN WEALTH POP PINTERG DENS INCOME ID GROWR;
run;

Data nym ; 
set sasuser.nym;
lex = log(expen);
lw = log(wealth);
lw2 = log(wealth)**2;
lpop = log(pop);
lpop2 = log(pop)**2;
lpop3 = log(pop)**3;
ldens = log(dens);
ldens2 = log(dens)**2;
lincome = log(income);
lincome2 = log(income)**2;

if growr < 0 then lgrowr = - log(1-0.2*growr);
if growr >= 0 then lgrowr =  log(1+0.2*growr);
IF pop < 32000;
IF pINTERG < 15;
IF GROWR > 0;
run;

proc reg data= nym;
model lEX = lW lW2 lPOP lPOP2 lPOP3 PINTERG lDENS lDENS2 lINCOME lINCOME2 lGROWR/ P R;
output out=b cookd=cd p=p r =r;
run;

proc gplot data =b;
plot cd*r R*p;
run;

DATA nym1; set nym; if _N_ ne 887 and _N_ ne 891;
;run;

proc reg data= nym;
model lEX = lW lW2 lPOP lPOP2 lPOP3 PINTERG lDENS lDENS2 lINCOME lINCOME2 lGROWR/ P R;
output out=b cookd=cd p=p r =r;
run;

proc gplot data =b;
plot cd*r R*p;
run;
proc reg data= nym;
model lEX = lW lW2 lPOP lPOP2 lPOP3 PINTERG lDENS lDENS2 lINCOME lINCOME2 lGROWR/ SELECTION=aDJRSQ;
run;

data boyle;
input V  P;
PM1 = 1/P;
PM2 = 1/P/P;
datalines;
    29.750  1.0
    19.125  1.5
    14.375  2.0
     9.500  3.0
     7.125  4.0
     5.625  5.0
     4.875  6.0
     4.250  7.0
     3.750  8.0
     3.375  9.0
     3.000 10.0
     2.625 12.0
     2.250 14.0
     2.000 16.0
     1.875 18.0
     1.750 20.0
     1.500 24.0
     1.375 28.0
     1.250 32.0
;


run;