\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlcom{#Monte Carlo study}
\hlkwd{set.seed}\hlstd{(}\hlnum{865177}\hlstd{)}
\hlkwd{library}\hlstd{(mvnfast)}
\hlstd{n}\hlkwb{<-}\hlnum{1000}
\hlstd{k}\hlkwb{<-}\hlnum{10}
\hlstd{Mu}\hlkwb{<-}\hlkwd{sample}\hlstd{(}\hlopt{-}\hlnum{5}\hlopt{:}\hlnum{5}\hlstd{,}\hlkwc{size}\hlstd{=k,}\hlkwc{rep}\hlstd{=T)}\hlopt{/}\hlnum{2}
\hlstd{Sigma}\hlkwb{<-}\hlnum{0.1}\hlopt{*}\hlnum{0.8}\hlopt{^}\hlkwd{abs}\hlstd{(}\hlkwd{outer}\hlstd{(}\hlnum{1}\hlopt{:}\hlstd{k,}\hlnum{1}\hlopt{:}\hlstd{k,}\hlstr{'-'}\hlstd{))}
\hlstd{al}\hlkwb{<-}\hlkwd{c}\hlstd{(}\hlnum{4}\hlstd{,}\hlkwd{rep}\hlstd{(}\hlnum{1}\hlstd{,}\hlnum{6}\hlstd{),}\hlkwd{rep}\hlstd{(}\hlopt{-}\hlnum{2}\hlstd{,}\hlnum{4}\hlstd{))}
\hlstd{be}\hlkwb{<-}\hlkwd{rep}\hlstd{(}\hlnum{0}\hlstd{,k}\hlopt{+}\hlnum{2}\hlstd{)}
\hlstd{be[}\hlnum{1}\hlopt{:}\hlnum{5}\hlstd{]}\hlkwb{<-}\hlkwd{c}\hlstd{(}\hlnum{2}\hlstd{,}\hlnum{1}\hlstd{,}\hlopt{-}\hlnum{2.0}\hlstd{,}\hlnum{2.2}\hlstd{,}\hlnum{3.6}\hlstd{)}
\hlstd{sigZ}\hlkwb{<-}\hlnum{5}
\hlstd{sigY}\hlkwb{<-}\hlnum{1}
\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{Z0}\hlkwb{<-}\hlkwd{round}\hlstd{(}\hlkwd{mean}\hlstd{(Z),}\hlnum{0}\hlstd{)}
\hlstd{zgrid}\hlkwb{<-}\hlkwd{seq}\hlstd{(Z0}\hlopt{-}\hlnum{5}\hlstd{,Z0}\hlopt{+}\hlnum{5}\hlstd{,}\hlkwc{by}\hlstd{=}\hlnum{0.1}\hlstd{)}
\hlstd{nZ}\hlkwb{<-}\hlkwd{length}\hlstd{(zgrid)}

\hlstd{nreps}\hlkwb{<-}\hlnum{200}
\hlstd{d}\hlkwb{<-}\hlnum{0.05}
\hlstd{ztrue.samp}\hlkwb{<-}\hlstd{zw.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}\hlopt{+}\hlnum{0.1}\hlopt{*}\hlstd{(Z}\hlopt{-}\hlstd{Z0)}\hlopt{^}\hlnum{2}\hlstd{,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{zest}\hlkwb{<-}\hlkwd{rep}\hlstd{(}\hlnum{0}\hlstd{,nZ)}
    \hlstd{ztrue}\hlkwb{<-}\hlstd{zw}\hlkwb{<-}\hlkwd{rep}\hlstd{(}\hlnum{0}\hlstd{,nZ)}
    \hlstd{fit.m}\hlkwb{<-}\hlstd{ps.fitted[,}\hlnum{1}\hlstd{]}
    \hlstd{fit.s}\hlkwb{<-}\hlstd{sigZ}
    \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,sigZ)}
        \hlstd{probu}\hlkwb{<-}\hlkwd{pnorm}\hlstd{(zu,fit.m,sigZ)}
        \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{ztrue[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{fitw}\hlkwb{<-}\hlkwd{lm}\hlstd{(Y}\hlopt{~}\hlstd{Z,}\hlkwc{subset}\hlstd{=(Z} \hlopt{>} \hlstd{zl} \hlopt{&} \hlstd{Z} \hlopt{<} \hlstd{zu),}\hlkwc{weight}\hlstd{=w)}
        \hlstd{zw[ix]}\hlkwb{<-}\hlkwd{sum}\hlstd{(}\hlkwd{coef}\hlstd{(fitw)} \hlopt{*} \hlkwd{c}\hlstd{(}\hlnum{1}\hlstd{,zgrid[ix]))}

    \hlstd{\}}
    \hlstd{ztrue.samp[irep,]}\hlkwb{<-}\hlstd{ztrue}
    \hlstd{zw.samp[irep,]}\hlkwb{<-}\hlstd{zw}
\hlstd{\}}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{+ Error in lm.wfit(x, y, w, offset = offset, singular.ok = singular.ok, : 0 (non-NA) cases}}\end{kframe}
\end{knitrout}