'%0 is the matrix of results
'%1 is the series,  
'%2 is the number of points (default=100)
'%3 is the left end of the distribution (minimum),
'%4 is the right end of the marginal distribution (maximum),
'       If %1=%2=0 then automatic selection/
'my_density_2001 price2001 

!pi=3.14159265358
     if %2="" then
          !M=100
     else
         !M=@val(%2)
     endif
matrix(!M,2) {%0}=0
!N=@obs({%1})            ' the number of observations
                                         '{%1} is the source series
!IQR=(@quantile({%1}, 0.75,"o")-@quantile({%1}, 0.25,"o"))/1.34
!s=@stdevp({%1})
     if !s<!IQR then            'Silverman bandwidth
         !h=0.9*!s/!N^0.2
     else
         !h=0.9*!IQR/!N^0.2
     endif
!min=@val(%3)
!max=@val(%4)
     if !min=0 or !min=NA then
          !min=@min({%1})-!h
     endif
     if !max=0 or !max=NA then
          !max=@max({%1})+!h
     endif

     for !i=1 to !M
        statusline i=!i
        {%0}(!i,1)=!min+(!i-1)*(!max-!min)/(!M-1)
          for !r=1 to !N
               {%0}(!i,2)={%0}(!i,2)+exp(-0.5*(({%0}(!i,1)-{%1}(!r))/!h)^2)/(!N*!h*sqr(2*!pi))
          next
     next
%0=%0+"_param"
table(10,4) {%0}
{%0}(1,1)="Series="
{%0}(1,2)=%1
{%0}(2,1)="min="
{%0}(2,2)=@min({%1})
{%0}(2,3)="max="
{%0}(2,4)=@max({%1})
{%0}(3,1)="s="
{%0}(3,2)=!s
{%0}(3,3)="IQR/1.34="
{%0}(3,4)=!IQR
{%0}(4,3)="IQR="
{%0}(4,4)=!IQR*1.34
{%0}(5,1)="h(s)="
{%0}(5,2)=0.9*!s/!N^0.2
{%0}(5,3)="h(IQR)="
{%0}(5,4)=0.9*!IQR/!N^0.2
{%0}(6,1)="h="
{%0}(6,2)=!h