cap program drop myspline
program myspline, sortpreserve
syntax varlist(max=1), Gen(string) [ Nknots(numlist) k(string)]
local v = "`varlist'"
if "`k'"=="" {
	if "`nknots'"=="" local nknots = 3
	sort `v'
	mat k = J(1,`nknots',.)
	forvalues i=1/`nknots' {
		mat k[1,`i'] = `v'[`i'/(`nknots'+1)*_N]
		}
	}
	else {
		mat k = `k'
		local nknots = colsof(k)
	}
gen `gen'1 = `v'
gen `gen'd1 = 1
forv i = 2/`=`nknots'-1' {
	/*
	local f = (k[1,`nknots']-k[1,`i'])/(k[1,`nknots']-k[1,1])
	gen `gen'`i' = max(0,`v'-k[1,`i'])^3-`f'*max(0,`v'-k[1,1])^3-(1-`f')*max(0,`v'-k[1,`nknots'])^3
	gen `gen'd`i' =  3*max(0,`v'-k[1,`i'])^2-`f'*3*max(0,`v'-k[1,1])^2-(1-`f')*3*max(0,`v'-k[1,`nknots'])^2
	*/
	// Stata code (see documentation mkspline)
	gen `gen'`i' = (max(0,`v'-k[1,`i'-1])^3- ///
		1/(k[1,`nknots']-k[1,`nknots'-1])*( ///
		max(0,`v'-k[1,`nknots'-1])^3*(k[1,`nknots']-k[1,`i'-1])- ///
		max(0,`v'-k[1,`nknots'])^3*(k[1,`nknots'-1]-k[1,`i'-1])))/ ///
		(k[1,`nknots']-k[1,1])^2
	gen `gen'd`i' = 3*(max(0,`v'-k[1,`i'-1])^2- ///
		1/(k[1,`nknots']-k[1,`nknots'-1])*( ///
		3*max(0,`v'-k[1,`nknots'-1])^2*(k[1,`nknots']-k[1,`i'-1])- ///
		3*max(0,`v'-k[1,`nknots'])^2*(k[1,`nknots'-1]-k[1,`i'-1])))/ ///
		(k[1,`nknots']-k[1,1])^2
	}
end

exit


clear
set obs 100
gen y = rnormal()
gen z = y^3 +rnormal()
mkspline t_ = y, nk(5) cubic
noi mat li r(knots)
mat a = r(knots)
myspline y, gen(s_) k(a)
noi mat li k
reg z s_?
predict p
reg z t_*
predict p1
sc z p p1 y
exit