\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{n}\hlkwb{<-}\hlnum{1000}
\hlstd{nreps}\hlkwb{<-}\hlnum{1000}
\hlstd{zest.samp}\hlkwb{<-}\hlkwd{matrix}\hlstd{(}\hlnum{0}\hlstd{,}\hlkwc{nrow}\hlstd{=nreps,}\hlkwc{ncol}\hlstd{=nZ)}
\hlstd{zest.dr.samp}\hlkwb{<-}\hlkwd{matrix}\hlstd{(}\hlnum{0}\hlstd{,}\hlkwc{nrow}\hlstd{=nreps,}\hlkwc{ncol}\hlstd{=nZ)}
\hlstd{zest.or.samp}\hlkwb{<-}\hlkwd{matrix}\hlstd{(}\hlnum{0}\hlstd{,}\hlkwc{nrow}\hlstd{=nreps,}\hlkwc{ncol}\hlstd{=nZ)}
\hlstd{gest.samp}\hlkwb{<-}\hlkwd{matrix}\hlstd{(}\hlnum{0}\hlstd{,}\hlkwc{nrow}\hlstd{=nreps,}\hlkwc{ncol}\hlstd{=nZ)}
\hlstd{gest2.samp}\hlkwb{<-}\hlkwd{matrix}\hlstd{(}\hlnum{0}\hlstd{,}\hlkwc{nrow}\hlstd{=nreps,}\hlkwc{ncol}\hlstd{=nZ)}
\hlstd{agest.samp}\hlkwb{<-}\hlkwd{matrix}\hlstd{(}\hlnum{0}\hlstd{,}\hlkwc{nrow}\hlstd{=nreps,}\hlkwc{ncol}\hlstd{=nZ)}
\hlkwa{for}\hlstd{(irep} \hlkwa{in} \hlnum{1}\hlopt{:}\hlstd{nreps)\{}
    \hlstd{X}\hlkwb{<-}\hlkwd{rmvn}\hlstd{(n,}\hlkwc{mu}\hlstd{=Mu,Sigma)}
    \hlstd{Xal}\hlkwb{<-}\hlkwd{cbind}\hlstd{(}\hlnum{1}\hlstd{,X)}
    \hlstd{ps.true}\hlkwb{<-}\hlstd{Xal} \hlopt{%*%} \hlstd{al}
    \hlstd{Z}\hlkwb{<-}\hlkwd{rnorm}\hlstd{(n,ps.true,sigZ)}
    \hlstd{Xm}\hlkwb{<-}\hlkwd{cbind}\hlstd{(}\hlnum{1}\hlstd{,Z,X)}
    \hlstd{muval}\hlkwb{<-}\hlkwd{as.vector}\hlstd{(Xm} \hlopt{%*%} \hlstd{be)}

    \hlstd{Y}\hlkwb{<-}\hlkwd{rnorm}\hlstd{(n,muval,sigY)}
    \hlstd{ps0}\hlkwb{<-}\hlkwd{lm}\hlstd{(Z}\hlopt{~}\hlstd{X)}
    \hlstd{ps.fitted}\hlkwb{<-}\hlkwd{predict}\hlstd{(ps0,}\hlkwc{newdata}\hlstd{=}\hlkwd{data.frame}\hlstd{(}\hlkwc{X}\hlstd{=X),}\hlkwc{interval}\hlstd{=}\hlstr{'prediction'}\hlstd{)}

    \hlstd{fit0}\hlkwb{<-}\hlkwd{lm}\hlstd{(Y}\hlopt{~}\hlstd{Z}\hlopt{+}\hlstd{X)}

    \hlstd{zest}\hlkwb{<-}\hlkwd{rep}\hlstd{(}\hlnum{0}\hlstd{,nZ)}
    \hlstd{zest.dr}\hlkwb{<-}\hlkwd{rep}\hlstd{(}\hlnum{0}\hlstd{,nZ)}
    \hlstd{zest.or}\hlkwb{<-}\hlkwd{rep}\hlstd{(}\hlnum{0}\hlstd{,nZ)}
    \hlstd{fit.m}\hlkwb{<-}\hlstd{ps.fitted[,}\hlnum{1}\hlstd{]}
    \hlstd{fit.s}\hlkwb{<-}\hlstd{(fit.m}\hlopt{-}\hlstd{ps.fitted[,}\hlnum{2}\hlstd{])}\hlopt{/}\hlnum{1.96}
    \hlstd{dvec}\hlkwb{<-}\hlkwd{dnorm}\hlstd{(Z,fit.m,fit.s)}

    \hlkwa{for}\hlstd{(ix} \hlkwa{in} \hlnum{1}\hlopt{:}\hlstd{nZ)\{}
        \hlstd{zl}\hlkwb{<-}\hlstd{zgrid[ix]}\hlopt{-}\hlstd{d}
        \hlstd{zu}\hlkwb{<-}\hlstd{zgrid[ix]}\hlopt{+}\hlstd{d}
        \hlstd{probl}\hlkwb{<-}\hlkwd{pnorm}\hlstd{(zl,fit.m,fit.s)}
        \hlstd{probu}\hlkwb{<-}\hlkwd{pnorm}\hlstd{(zu,fit.m,fit.s)}
        \hlstd{probix}\hlkwb{<-}\hlstd{probu}\hlopt{-}\hlstd{probl}
        \hlstd{w}\hlkwb{<-}\hlstd{(Z} \hlopt{>} \hlstd{zl} \hlopt{&} \hlstd{Z} \hlopt{<} \hlstd{zu)}\hlopt{/}\hlstd{probix}
        \hlstd{w}\hlkwb{<-}\hlstd{w}\hlopt{*}\hlstd{Z}\hlopt{/}\hlstd{dvec}
        \hlstd{zest[ix]}\hlkwb{<-}\hlkwd{sum}\hlstd{(w}\hlopt{*}\hlstd{Y,}\hlkwc{na.rm}\hlstd{=T)}\hlopt{/}\hlkwd{sum}\hlstd{(w,}\hlkwc{na.rm}\hlstd{=T)}
        \hlstd{fit1}\hlkwb{<-}\hlkwd{lm}\hlstd{(Y}\hlopt{~}\hlstd{Z)}
        \hlstd{mu.fit}\hlkwb{<-}\hlkwd{predict}\hlstd{(fit1,}\hlkwd{data.frame}\hlstd{(}\hlkwc{Z}\hlstd{=}\hlkwd{rep}\hlstd{(zgrid[ix],n)))}
        \hlstd{zest.dr[ix]}\hlkwb{<-}\hlkwd{sum}\hlstd{(w}\hlopt{*}\hlstd{(Y}\hlopt{-}\hlstd{mu.fit),}\hlkwc{na.rm}\hlstd{=T)}\hlopt{/}\hlkwd{sum}\hlstd{(w,}\hlkwc{na.rm}\hlstd{=T)}\hlopt{+}\hlkwd{mean}\hlstd{(mu.fit)}
        \hlstd{zest.or[ix]}\hlkwb{<-}\hlkwd{mean}\hlstd{(mu.fit)}
    \hlstd{\}}
    \hlstd{zest.samp[irep,]}\hlkwb{<-}\hlstd{zest}
    \hlstd{zest.dr.samp[irep,]}\hlkwb{<-}\hlstd{zest.dr}
    \hlstd{zest.or.samp[irep,]}\hlkwb{<-}\hlstd{zest.or}
    \hlstd{aw}\hlkwb{<-}\hlkwd{abs}\hlstd{(Z}\hlopt{-}\hlstd{fit.m)}
    \hlstd{afit}\hlkwb{<-}\hlkwd{lm}\hlstd{(Y}\hlopt{~}\hlstd{Z,} \hlkwc{weight}\hlstd{=aw)}
    \hlstd{agest.samp[irep,]}\hlkwb{<-}\hlkwd{cbind}\hlstd{(}\hlnum{1}\hlstd{,zgrid)} \hlopt{%*%} \hlkwd{coef}\hlstd{(afit)}

    \hlstd{gfit}\hlkwb{<-}\hlkwd{lm}\hlstd{(Y}\hlopt{~}\hlstd{Z}\hlopt{+}\hlstd{fit.m)}
    \hlstd{gest.samp[irep,]}\hlkwb{<-}\hlkwd{cbind}\hlstd{(}\hlnum{1}\hlstd{,zgrid,}\hlkwd{mean}\hlstd{(fit.m))} \hlopt{%*%} \hlkwd{coef}\hlstd{(gfit)}


    \hlstd{Z1}\hlkwb{<-}\hlkwd{cbind}\hlstd{(}\hlnum{1}\hlstd{,Z}\hlopt{-}\hlstd{fit.m)}
    \hlstd{Z2}\hlkwb{<-}\hlkwd{cbind}\hlstd{(}\hlnum{1}\hlstd{,Z)}
    \hlstd{g2fit.est}\hlkwb{<-}\hlkwd{solve}\hlstd{(}\hlkwd{t}\hlstd{(Z1)} \hlopt{%*%} \hlstd{Z2)} \hlopt{%*%} \hlstd{(}\hlkwd{t}\hlstd{(Z1)} \hlopt{%*%} \hlstd{Y)}
    \hlstd{gest2.samp[irep,]}\hlkwb{<-}\hlkwd{cbind}\hlstd{(}\hlnum{1}\hlstd{,zgrid)} \hlopt{%*%} \hlstd{g2fit.est}


\hlstd{\}}
\end{alltt}
\end{kframe}
\end{knitrout}