#delim ;
version 6.0;
*
 Certify bspline and frencurv using the Auto data
*;

* Open log file *;log using cert1.log,replace;
capture noisily{
* Beginning of main capture block *;

* Begin certification script *;
cscript "Certification of bspline and frencurv"
 adofiles bspline frencurv;

* Use auto data set and define additional variables *;
use auto,clear;sort weight;desc;

* Summarize weight (which will be the X-variate) *;
summ weight,detail;return list;

* Error conditions for bspline *;
* Negative power *;
rcof "noisily bspline,x(weight) k(1760(770)4840) gen(bs) pow(-2)"==498;
* No observations *;
rcof "noisily bspline if(0),x(weight) k(1760(770)4840) gen(bs) pow(3)"==2000;
* No varlist or generate prefix *;
rcof "noisily bspline,x(weight) k(1760(770)4840) pow(3)"==498;
* Generate prefix too long for sensible spline names *;
rcof "noisily bspline,x(weight) k(1760(770)4840) gen(abcdefgh) pow(3)"==198;
* Generate prefix too long for sensible knot names *;
rcof "noisily bspline,x(weight) k(1760(770)4840) gen(abcdefg) pow(3)"==198;
* Too many knots for default matsize *;
rcof "noisily bspline,x(weight) k(1760(28)4840) gen(bs) pow(3)"==908;
* Too many knots for Stata to handle in a matrix define expression *;
rcof "noisily bspline,x(weight) k(1760(7)4840) gen(bs) pow(3)"==130;
* Show that the standard version of the command works *;
noisily bspline,x(weight) k(1760(770)4840) gen(bs) pow(3);
regr mpg bs*,noconst robust;
drop bs*;
return list;

* Error conditions for frencurv *;
* Negative power *;
rcof "noisily frencurv,x(weight) r(1760(770)4840) gen(sp) pow(-2)"==498;
* No observations *;
rcof "noisily frencurv if(0),x(weight) r(1760(770)4840) gen(sp) pow(3)"==2000;
* No varlist or generate prefix *;
rcof "noisily frencurv,x(weight) r(1760(770)4840) pow(3)"==498;
* Generate prefix too long for sensible spline names *;
rcof "noisily frencurv,x(weight) r(1760(770)4840) gen(abcdefgh) pow(3)"==198;
* Too many reference values for default matsize *;
rcof "noisily frencurv,x(weight) r(1760(28)4840) gen(sp) pow(3)"==908;
* Too many reference values for Stata to handle in a matrix define expression *;
rcof "noisily frencurv,x(weight) r(1760(7)4840) gen(sp) pow(3)"==130;
* Reference points out of range of knots *;
rcof "noisily frencurv,x(weight) r(1760(770)4840) k(1760(770)4840) noexk gen(sp) pow(3)"==498;
* Reference points do not identify B-splines *;
rcof "noisily frencurv,x(weight) r(1760(1)1764) k(1760(770)4840) gen(sp) pow(3)"==498;
* Show that the standard version of the command works *;
frencurv,x(weight) r(1760(770)4840) gen(sp) pow(3);
return list;desc sp*;
regr mpg sp*,noconst robust;
drop sp*;

*
 Show that in and if work as expected
 (alone and in combination)
 in that generated splines have as many valid values as the x-variate
 on which they were based
*;
summ weight if(!foreign);local nweius=r(N);
summ weight in 1/37;local nwei37=r(N);
summ weight if(!foreign) in 25/l;local nweiusl=r(N);
* Test bspline *;
bspline if(!foreign),x(weight) k(1760(770)4840) gen(bsp) pow(1);
summ bsp*;assert r(N)==`nweius';
drop bsp*;
bspline in 1/37,x(weight) k(1760(770)4840) gen(bsp) pow(1);
summ bsp*;assert r(N)==`nwei37';
drop bsp*;
bspline if(!foreign) in 25/l,x(weight) k(1760(770)4840) gen(bsp) pow(1);
summ bsp*;assert r(N)==`nweiusl';
drop bsp*;
* Test frencurv *;
frencurv if(!foreign),x(weight) r(1760(770)4840) gen(sp) pow(1);
summ sp*;assert r(N)==`nweius';
drop sp*;
frencurv in 1/37,x(weight) r(1760(770)4840) gen(sp) pow(1);
summ sp*;assert r(N)==`nwei37';
drop sp*;
frencurv if(!foreign) in 25/l,x(weight) r(1760(770)4840) gen(sp) pow(1);
summ sp*;assert r(N)==`nweiusl';
drop sp*;

*
 Show that frencurv still works with even powers
 and knots off-centre with respect to reference value midpoints
*;
frencurv,x(weight) r(1760(770)4840) k(1560(770)5410) gen(sp) pow(0);
desc sp*;
regr mpg sp*,noconst robust;
drop sp*;
frencurv,x(weight) r(1760(770)4840) k(1560(770)5410) gen(sp) pow(2);
desc sp*;
regr mpg sp*,noconst robust;
drop sp*;
frencurv,x(weight) r(1760(770)4840) k(1560(770)5410) gen(sp) pow(4);
desc sp*;
regr mpg sp*,noconst robust;
drop sp*;

*
 Show that bspline and frencurv give the same result as mkspline
 if power of spline is 1
*;
scal tiny=1e-6;
mkspline msp1 1760 msp2 2530 msp3 3300 msp4 4070 msp5=weight;
regr mpg msp*,noconst robust;
predict mpghatm;
frencurv,x(weight) r(1760,2530,3300,4070,4840) p(1) g(rsp);
regr mpg rsp*,noconst robust;
predict mpghatr;
bspline,x(weight) k(1760,2530,3300,4070,4840) p(1) g(bsp);
regr mpg bsp*,noconst robust;
predict mpghatb;
assert reldif(mpghatm,mpghatb)<tiny;
assert reldif(mpghatm,mpghatr)<tiny;
assert reldif(mpghatb,mpghatr)<tiny;
drop mpghat* msp* bsp* rsp*;

*
 Show that bspline and frencurv give the same result
 as the plus-function basis
 for each choice of positive power up to 3
 (after which I find that there are stability problems)
*;
scal tiny=1e-6;
local rlist "1760(770)4840";
local powmax=3;local power=-1;
while(`power'<`powmax'){local power=`power'+1;
  frencurv,x(weight) r(`rlist') g(rsp) type(double)
    pow(`power');
  local knots "`r(knots)'";local nknot:word count `knots';
  bspline,x(weight) k(`knots') noexk g(bsp) type(double)
    pow(`power');
  local i1=0;
  while(`i1'<`nknot'-`power'-1){local i1=`i1'+1;
    local origin:word `i1' of `knots';
    gene double pf`i1'=(weight-`origin')^`power';
    replace pf`i1'=0 if(weight<`origin');
    label variable pf`i1' "Plus-function at `origin'";
  };
  disp "Reference splines of power `power'";
  desc rsp*;regr mpg rsp*,noconst robust;predict mpghatr;
  disp "B-splines of power `power'";
  desc bsp*;regr mpg bsp*,noconst robust;predict mpghatb;
  disp "Plus-functions of power `power'";
  desc pf*;regr mpg pf*,noconst robust;predict mpghatp;
  disp "Reference spline vs B-spline fitted values of power `power'";
  * Check that fitted values are the same for reference splines and B-splines *;
  plot mpghatr mpghatb;
  assert reldif(mpghatr,mpghatb)<tiny;
  *
   Check that fitted values are the same
   for reference splines and plus-functions
  *;
  disp "Reference spline vs plus-function spline fitted values of power `power'";
  plot mpghatr mpghatp;
  assert reldif(mpghatr,mpghatp)<tiny;
  * Check that fitted values are the same for B-splines and plus-functions *;
  disp "B-spline vs plus-function spline fitted values of power `power'";
  plot mpghatb mpghatp;
  assert reldif(mpghatb,mpghatp)<tiny;
  drop mpghat* rsp* bsp* pf*;
};

*
 Show that bspline and frencurv give the same result as each other
 for each choice of positive power up to 7
 (after which I find that there are stability problems)
*;
scal tiny=1e-6;
local rlist "1760(770)4840";
local powmax=7;local power=-1;
while(`power'<`powmax'){local power=`power'+1;
  frencurv,x(weight) r(`rlist') g(rsp) type(double)
    pow(`power');
  local knots "`r(knots)'";local nknot:word count `knots';
  bspline,x(weight) k(`knots') noexk g(bsp) type(double)
    pow(`power');
  disp "Reference splines of power `power'";
  desc rsp*;regr mpg rsp*,noconst robust;predict mpghatr;
  disp "B-splines of power `power'";
  desc bsp*;regr mpg bsp*,noconst robust;predict mpghatb;
  disp "Reference spline vs B-spline fitted values of power `power'";
  plot mpghatr mpghatb;
  assert reldif(mpghatr,mpghatb)<tiny;
  drop mpghat* rsp* bsp*;
};


*
 Show that bspline and frencurv give the same result as each other
 if reference points are all within realized range of X-variate
 for each choice of positive power up to 10
 (after which there are stability problems)
*;
scal tiny=1e-6;
local klist "1760(770)4840";
local powmax=10;local power=0;
while(`power'<`powmax'){local power=`power'+1;
  local space=(4840-1760)/(4+`power'-1);
  local rlist "1760(`space')4840";
  frencurv,x(weight) r(`rlist') noexr k(`klist') g(rsp) type(double)
    pow(`power');
  bspline,x(weight) k(`klist') g(bsp) type(double)
    pow(`power');
  disp "Reference splines of power `power'";
  desc rsp*;regr mpg rsp*,noconst robust;predict mpghatr;
  *
   Show that maximum and minimum weights correspond to fitted values
   equal to top and bottom entries of the estimate vector, respectively
   (recall that data set is sorted by weight,
   and that minimum and maximum reference points correspond
   to minimum and maximum weights)
  *;
  matr def beta=e(b);local nparm=colsof(beta);
  scal spmin=beta[1,1];scal spmax=beta[1,`nparm'];
  assert reldif(spmin,mpghatr)<tiny if(_n==1);
  assert reldif(spmax,mpghatr)<tiny if(_n==_N);
  * Test that fitted values are the same for B-spline and reference spline *;
  disp "B-splines of power `power'";
  desc bsp*;regr mpg bsp*,noconst robust;predict mpghatb;
  disp "Reference spline vs B-spline fitted values of power `power'";
  plot mpghatr mpghatb;
  assert reldif(mpghatr,mpghatb)<tiny;
  drop mpghat* rsp* bsp*;
};

* End of main capture block *;
};
* Close log file and exit *;log close;exit;