/* ================================================================================== */ /* Procedure for: Locally Linear Estimation with Plug-In Bandwidth */ /* Stefan Hupfeld, August 2006 */ /* The Plug-In Method follows Loader (1999, 2004) */ /* ================================================================================== */ /* http://creativecommons.org/licenses/by-nc-sa/2.0/de/ */ /* Under the above link you find the license for this program. */ /* You are free: -to copy, distribute, display, and perform the work */ /* -to make derivative works. */ /* Attribution. You must attribute the work in the manner specified by the author or licensor. */ /* Noncommercial. You may not use this work for commercial purposes. */ /* Share Alike. If you alter, transform, or build upon this work, */ /* you may distribute the resulting work only under a license identical to this one. */ /* For any reuse or distribution, you must make clear to others the license terms of this work. */ /* Any of these conditions can be waived if you get permission from the copyright holder. */ /* =====HEAD========================================================================= */ new; library pgraph; cls; graphset; pqgwin many; /* =====MAIN========================================================================= */ n = 100; e = rndn(n,1); x = seqa(1,1,n); y = -0.1*(x-1)^2+0.01*x^3+500*e; /* ==PILOT=========================================================================== */ bw = local_quad("Pilot Estimation", x, y); /* ==LL-Estimation==================================================================== */ y_fit = zeros(n, 1); y_fit = local_lin("Local Linear Estimation", bw, x, y); plot2d("Local Linear Estimation", "x", "y, y_fit", x, y_fit~y); ndpclex; end; /* ================================================================================== */ /* ========================================================================================== */ /* Procedures */ /* ========================================================================================== */ /* ================================================================================== */ /* 2D-Plot */ /* Input: String for Title, x-Axis, y-Axis, and xy-Values */ /* Output: Standardized xy-Plot */ /* Stefan Hupfeld, Oct. 14 2005 */ /* ================================================================================== */ proc plot2d(plot_title, plot_x_label, plot_y_label, plot_x, plot_y); _pcolor={0 7 0 7}; _pmcolor = {0, 0, 0, 0, 0, 0, 0, 0, 15}; _pltype = { 0, 0, 0, 0 }; _paxht = 0.2; _ptitlht = 0.2; title(plot_title); xlabel(plot_x_label); ylabel(plot_y_label); xy(plot_x, plot_y); retp("(Plot)"); endp; /* ================================================================================== */ /* Local Linear Estimation */ /* Input: Estimation Title, Bandwidth, x-Values, y-Values */ /* Output: LL-Fit */ /* Stefan Hupfeld, August 2006 */ /* ================================================================================== */ proc local_lin(est_title, bw, x_val, y_val); local y_ll_fit, nn, x0, beta, w, XWLS; nn = rows(x_val); y_ll_fit = zeros(nn, 1); output off; format /mat /ld 6,0; for i (1, nn, 1); locate(1,1); print /flush est_title " " i "of " nn; x0 = x_val[i,1]; w = knorm(x_val, x0, bw); XWLS = ones(nn, 1)~(x_val - x0); beta = inv((XWLS.*w)'*XWLS)*(XWLS.*w)'*y_val; y_ll_fit[i] = beta[1]; endfor; format /mat /ld 6,3; output on; retp(y_ll_fit); endp; /* ================================================================================== */ /* Local Quad Estimation */ /* Input: Estimation Title, x-Values, y-Values */ /* Output: bw_star */ /* Stefan Hupfeld, August 2006 */ /* ================================================================================== */ proc local_quad(est_title, x_val, y_val); local nq, bwq, yq, x0q, H1q, H2q, e_q, wq, XWLSq, betaq, holdq, bw_star; nq = rows(x_val); bwq = 5*nq^(-1/5); @ Pilot BW; Order of Magn ~n^(-1/5) @ print "bw pilot" bwq; yq = zeros(nq, 3); @ 3 for 2nd Derivative @ e_q = zeros(1,3); @ 3 for 2nd Derivative @ e_q[1] = 1; H1q = zeros(nq, 1); H2q = zeros(nq, 1); output off; format /mat /ld 6,0; for i (1, nq, 1); locate(1,1); print /flush est_title i " of " nq; x0q = x_val[i,1]; wq = knorm(x_val, x0q, bwq); XWLSq = ones(nq, 1)~(x_val - x0q)~(x_val - x0q)^2; betaq = inv((XWLSq.*wq)'*XWLSq)*(XWLSq.*wq)'*y_val; yq[i,.] = betaq'; holdq = e_q*inv((XWLSq.*wq)'*XWLSq)*(XWLSq.*wq)'; H1q[i] = holdq[i]; H2q[i] = sqrt(sumc((holdq^2)')); endfor; format /mat /ld 6,3; output on; bw_star = bwidth(H1q, H2q, yq, y_val, x_val); print "bw_star " est_title bw_star; retp(bw_star); endp; /* ================================================================================== */ /* Plug-In Bandwidth Loader 2004 */ /* Source: Own, Jan 2006 */ /* ================================================================================== */ proc bwidth(H1_b, H2_b, y_pilot_b, dep_b, x_b); local nu1_b, nu2_b, sigma_b, K1_b, K2_b, frac1_b, frac2_b, sum_der_b, n_b, bwidth; n_b = rows(x_b); K1_b = 0.282095; K2_b = 1; nu1_b = sumc(H1_b); nu2_b = sumc(H2_b^2); sigma_b = sqrt(sumc((dep_b-y_pilot_b[.,1])^2)*(n_b-2*nu1_b+nu2_b)^(-1)); y_pilot_b[.,3] = y_pilot_b[.,3]*2; sum_der_b = sumc(y_pilot_b[.,3]^2); frac1_b = sigma_b*(maxc(x_b)-minc(x_b))*K1_b; frac2_b = n_b*K2_b*sum_der_b; bwidth = (frac1_b/frac2_b)^(-1/5); retp(bwidth); endp; /* ================================================================================== */ /* Normal Kernel */ /* ================================================================================== */ proc knorm(xx, x0, h); local k,n; k = cols(xx); n = rows(xx); retp(exp(-(0.5/h^2)*sumc(((xx-x0).*(xx-x0))'))/(sqrt(2*pi*h^2)^k)); endp;