data one;
input x y;
cards;
1950	604.892513
1951	638.4976526
1952	1243.390166
1953	806.5233506
1954	1008.154188
1955	537.6822338
1956	1108.969607
1957	840.1284902
1958	1747.46726
1959	1344.205584
1960	1747.46726
1961	1613.046701
1962	1276.995305
1963	1411.415864
1964	1545.836422
1965	1881.887818
1966	1411.415864
1967	1142.574747
1968	1545.836422
1969	1881.887818
1970	1747.46726
1971	1344.205584
1972	1545.836422
1973	2016.308377
1974	1411.415864
1975	1613.046701
1976	1613.046701
1977	1814.677539
1978	1814.677539
1979	2553.99061
1980	2016.308377
1981	1814.677539
1982	2217.939214
1983	2352.359773
1984	2419.570052
1985	2016.308377
1986	1949.098097
1987	1814.677539
1988	2419.570052
1989	1814.677539
1990	2217.939214
1991	1814.677539
1992	1949.098097
1993	2016.308377
1994	1814.677539
1995	1411.415864
1996	1276.995305
1997	2150.728935
1998	2553.99061
1999	2352.359773
2000	2285.149493
2001	2217.939214
2002	1881.887818
2003	2621.20089
2004	2352.359773
2005	2150.728935
2006	1613.046701
2007	1881.887818
2008	2486.780331
2009	1478.626143
2010	2083.518656
2011	1478.626143
2012	2419.570052
;
proc nlin data = one best = 3;
  parms b0=500 to 700 by 20 b1=1950 to 2010 by 2 njoint=1975 to 1985 by 1;
  if x<njoint then do;
  model y = b0 + b1*x;
der.b0=1;
der.b1=x;
der.njoint=0;
end;
else do;
model y=b0+b1*njoint;
der.b0=1;
der.b1=njoint;
der.njoint=b1;
end;
file print;
if _obs_ =1 and _model_ =0 then do;
  plateau = b0 + b1*njoint;
  put plateau=;
end;
  plateau=b0+b1*njoint;
id plateau;
  output out = new p = pry parms=b0 b1 njoint sse=sse;
run;
proc plot;
plot y*x='+' pry*x='*'/overlay;
run;
proc means noprint;
var y sse b0 b1 njoint plateau;
output out = new2 n = tdf
mean = y sse b0 b1 njoint plateau
css=csst;
data new3; set new2;
intercpt=b0; slope=b1; joint=njoint;
rsq=(csst-sse)/csst;
edf=tdf-3;
ssr=csst-sse;
msr=ssr/2;
mse=sse/edf;
f=msr/mse;
probf=1-(probf(f,2,edf));
keep intercpt slope joint plateau rsq f probf;
proc print;
run;