#delim ;
version 17.0;
*
 Fit a ridit spline model of npm with respect to  tons
 in the US cars of the xauto dataset as the training set,
 and use the model to predict npm
 in the non-US cars of the auto dataset as the test set.
 SSC packages required:
 polyspline, parmest, wridit, eclplot
*;

clear all;
set scheme rbn4mono;


*
 Input US cars in the xauto data as training set
*;
qui xauto, clear;
keep if us==1;
describe, full;


*
 Define native ridit of tons
*;
wridit tons, gene(r_tons);
lab var r_tons "Ridit of tons in US cars";
summ r_tons, de format;


*
 Define ridit spline basis in tons
*;
polyspline r_tons, power(3) refpts(0(0.25)1)
  gene(sp_) labprefix("Spline at ");
describe sp_*, fu;


*
 Fit ridit spline model
*;
regress npm sp_*, vce(robust) noconst;
parmest, label escal(N) rename(es_1 N) format(estimate min* max* %8.3f)
  list(N parm label estimate min* max*, abbr(32));


*
 Display real and predicted values of npm in training set
*;
predict npmhat, xb;
predict senpmhat, stdp;
gene dof=e(df_r);
lab var npmhat "Predicted fuel consumption";
lab var senpmhat "SE of predicted fuel consumption";
lab var dof "Degrees of freedom";
parmcip, estimate(npmhat) stderr(senpmhat);
describe, full;
eclplot npmhat min* max* tons, eplot(line) rplot(rline)
  estopts(sort) ciopts(sort)
  addplot(scatter npm tons, msym(circle) msize(1))
  ylabel(5(1)25)
  ytitle("Fuel consumption (nipperkins/mile) in training set of US cars")
  xlabel(0.75(0.25)2.5)
  xtitle("Weight (US tons)")
  legend(on size(3) row(1) order(3 "Observed" 1 "Predicted (with 95% CI)"))
  ysize(4.6) xsize(4.6);
graph export figseq3.pdf, replace;
more;


*
 Save estimation results and frequencies of tons
*;
estimates save rstons.ster, replace;
xcontract tons, saving(tonsfreq.dta, replace)
  list(, abbr(32));


*
 Input non-US cars in the xauto data as testg set
*;
qui xauto, clear;
keep if us==0;
describe, full;


*
 Define foreign ridit of tons
*;
cap frame drop tonsfreq;
frame create tonsfreq;
frame tonsfreq: use tonsfreq.dta, clear;
fridit tons, fframe(tonsfreq, weight(_freq))
  gene(r_tons);
lab var r_tons "Ridit of tons in US cars";
summ r_tons, de format;


*
 Define ridit spline basis in tons
*;
polyspline r_tons, power(3) refpts(0(0.25)1)
  gene(sp_) labprefix("Spline at ");
describe sp_*, fu;


*
 Use ridit spline model from training set
*;
estimates use rstons.ster;
regress;
parmest, label escal(N) rename(es_1 N) format(estimate min* max* %8.3f)
  list(N parm label estimate min* max*, abbr(32));


*
 Display real and predicted values of npm in test set
*;
predict npmhat, xb;
predict senpmhat, stdp;
gene dof=e(df_r);
lab var npmhat "Predicted fuel consumption";
lab var senpmhat "SE of predicted fuel consumption";
lab var dof "Degrees of freedom";
parmcip, estimate(npmhat) stderr(senpmhat);
describe, full;
eclplot npmhat min* max* tons, eplot(line) rplot(rline)
  estopts(sort) ciopts(sort)
  addplot(scatter npm tons, msym(circle) msize(1))
  ylabel(5(1)25)
  ytitle("Fuel consumption (nipperkins/mile) in test set of non-US cars")
  xlabel(0.75(0.25)2.5)
  xtitle("Weight (US tons)")
  legend(on size(3) row(1) order(3 "Observed" 1 "Predicted (with 95% CI)"))
  ysize(4.6) xsize(4.6);
graph export figseq4.pdf, replace;
more;


exit;