/* Program to accompany the articles "Smoothers for periodic data" (http://blogs.sas.com/content/iml/2012/04/05/smoothers-for-periodic-data/) and "Creating a periodic smoother" (http://blogs.sas.com/content/iml/2012/04/06/creating-a-periodic-smoother/) by Rick Wicklin Create a periodic loess smoother. */ /* Read the data directly from the Internet */ filename webfile url "http://academic.udayton.edu/kissock/http/Weather/gsod95-current/NYALBANY.txt" /* proxy='http://...:80' */; data TempData; infile webfile; input month day year Temperature; format date date9.; date=MDY(month, day, year); dayofyear=0; dayofyear=put(date,julday3.); /* incorporate leap years into calculations */ Proportion = dayofyear / put(mdy(12, 31, year), julday3.); label Proportion = "Day of Year (Normalized)"; CurrentWinter = (date >='01dec2011'd & date<='15mar2012'd); if Temperature^=-99 then output; run; data TempData; set TempData TempData(where=(CurrentWinter=1) rename=(Proportion=Prop Temperature=Temp)); run; proc sgplot data=TempData; scatter x=Proportion y=Temperature / transparency=0.8; scatter x=Prop y=Temp / markerattrs=(color='gray' symbol=CircleFilled) legendlabel="Winter 2011-12"; loess x=Proportion y=Temperature / smooth=0.167 nomarkers lineattrs=(thickness=4px) legendlabel="Loess smoother"; yaxis grid; title "Temperature in Albany, NY (1995-2012)"; run; /* where did smoothing parameter come from? */ /* smooth = 61/365 = 0.167 ==> each day predicted by local weighted regression with +-30 days */ /* proc loess data=TempData plots(only maxpoints=none)=(FitPlot CriterionPlot); model Temperature = Proportion/ smooth=0.167 interp=cubic; run; */ /* extend data to each side */ data Periodic; set TempData(in=before) TempData TempData(in=after); if before then proportion = proportion - 1; /* (-1,0] */ if after then proportion = proportion + 1; /* (1,2] */ run; /* 3 times the data; use 1/3 of smoothing parameter */ data Score; do proportion = 0 to 1 by 1/365; output; end; proc loess data=Periodic plots(only maxpoints=none)=(FitPlot CriterionPlot); model Temperature = Proportion/ smooth=0.0557 interp=cubic; score data=Score; ods output ScoreResults=Fit; run; data Combine; merge TempData Fit(rename=(Proportion=Prop)); run; proc sgplot data=Combine; scatter x=Proportion y=Temperature / transparency=0.8; scatter x=Prop y=Temp / markerattrs=(color='gray' symbol=CircleFilled) legendlabel="Winter 2011-12"; series x=Prop y=P_Temperature/ lineattrs=(thickness=4px) legendlabel="Periodic smoother"; /* truly periodic */ yaxis grid; title "Temperature in Albany, NY (1995-2012)"; run;