## "Fossies" - the Fresh Open Source Software Archive

### Member "gretl-2019c/doc/tex/gretlR.tex" (8 Jan 2019, 22335 Bytes) of package /linux/misc/gretl-2019c.tar.xz:

As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) TeX and LaTeX source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file.

    1 \chapter{Gretl and R}
2 \label{chap:gretlR}
3
4 \section{Introduction}
5 \label{R-intro}
6
7 \app{R} is, by far, the largest free statistical
8 project.\footnote{\app{R}'s homepage is at
9   \url{http://www.r-project.org/}.} Like gretl, it is a GNU project
10 and the two have a lot in common; however, gretl's approach focuses on
11 ease of use much more than \app{R}, which instead aims to encompass
12 the widest possible range of statistical procedures.
13
14 As is natural in the free software ecosystem, we don't view ourselves
15 as competitors to \app{R},\footnote{OK, who are we kidding? But it's
16   \emph{friendly} competition!} but rather as projects sharing a
17 common goal who should support each other whenever possible. For this
18 reason, gretl provides a way to interact with \app{R} and thus enable
19 users to pool the capabilities of the two packages.
20
21 In this chapter, we will explain how to exploit \app{R}'s power from
22 within gretl. We assume that the reader has a working installation of
23 \app{R} available and a basic grasp of \app{R}'s syntax.\footnote{The
24   main reference for \app{R} documentation is
26   tutorials abound on the Net; as always, Google is your friend.}
27
28 Despite several valiant attempts, no graphical shell has gained wide
29 acceptance in the \app{R} community: by and large, the standard method
30 of working with \app{R} is by writing scripts, or by typing commands
31 at the \app{R} prompt, much in the same way as one would write
32 gretl scripts or work with the gretl console. In this
33 chapter, the focus will be on the methods available to execute \app{R}
34 commands without leaving gretl.
35
36 \section{Starting an interactive \app{R} session}
37 \label{sec:R-interactive}
38
39 The easiest way to use \app{R} from gretl is in interactive
40 mode.  Once you have your data loaded in gretl, you can select
41 the menu item Tools, Start GNU R'' and an interactive \app{R}
43
44 \subsection{A simple example: OLS on cross-section data}
45 \label{sec:R-ols-ex}
46
47 For this example we use Ramanathan's dataset \texttt{data4-1}, one of
48 the sample files supplied with gretl.  We first run, in
49 gretl, an OLS regression of \texttt{price} on \texttt{sqft},
50 \texttt{bedrms} and \texttt{baths}.  The basic results are shown in
51 Table \ref{tab:data4-1-gretlOLS}.
52
53 \begin{table}[htbp]
54 \caption{OLS house price regression via gretl}
55 \label{tab:data4-1-gretlOLS}
56 \begin{center}
57
58 \begin{tabular*}{0.75\textwidth}{@{\extracolsep{\fill}}
59 l% col 1: varname
60   D{.}{.}{-1}% col 2: coeff
61     D{.}{.}{-1}% col 3: sderr
62       D{.}{.}{-1}% col 4: t-stat
63         D{.}{.}{4}}% col 5: p-value (or slope)
64 Variable &
65   \multicolumn{1}{c}{Coefficient} &
66     \multicolumn{1}{c}{Std.\ Error} &
67       \multicolumn{1}{c}{$t$-statistic} &
68         \multicolumn{1}{c}{p-value} \$1ex] 69 const & 70 129.062 & 71 88.3033 & 72 1.4616 & 73 0.1746 \\ 74 sqft & 75 0.154800 & 76 0.0319404 & 77 4.8465 & 78 0.0007 \\ 79 bedrms & 80 -21.587 & 81 27.0293 & 82 -0.7987 & 83 0.4430 \\ 84 baths & 85 -12.192 & 86 43.2500 & 87 -0.2819 & 88 0.7838 \\ 89 \end{tabular*} 90 \end{center} 91 \end{table} 92 93 We will now replicate the above results using \app{R}. Select 94 the menu item Tools, Start GNU R''. A window similar to the one 95 shown in figure \ref{fig:Rwind1} should appear. 96 97 \begin{figure}[htbp] 98 \centering 99 \includegraphics[scale=0.7]{figures/Rwindow-1} 100 \caption{\app{R} window} 101 \label{fig:Rwind1} 102 \end{figure} 103 104 The actual look of the \app{R} window may be somewhat different from 105 what you see in Figure~\ref{fig:Rwind1} (especially for Windows 106 users), but this is immaterial. The important point is that you have a 107 window where you can type commands to \app{R}. If the above procedure 108 doesn't work and no \app{R} window opens, it means that gretl 109 was unable to launch \app{R}. You should ensure that \app{R} is 110 installed and working on your system and that gretl knows where 111 it is. The relevant settings can be found by selecting the Tools, 112 Preferences, General'' menu entry, under the Programs'' tab. 113 114 Assuming \app{R} was launched successfully, you will see notification 115 that the data from gretl are available. In the background, 116 gretl has arranged for two \app{R} commands to be executed, one 117 to load the gretl dataset in the form of a \emph{data frame} 118 (one of several forms in which \app{R} can store data) and one to 119 \emph{attach} the data so that the variable names defined in the 120 gretl workspace are available as valid identifiers within 121 \app{R}. 122 123 In order to replicate gretl's OLS estimation, go into the \app{R} 124 window and type at the prompt 125 \begin{code} 126 model <- lm(price ~ sqft + bedrms + baths) 127 summary(model) 128 \end{code} 129 130 \begin{figure}[htbp] 131 \centering 132 \includegraphics[scale=0.7]{figures/Rwindow-2} 133 \caption{OLS regression on house prices via \app{R}} 134 \label{fig:Rwind2} 135 \end{figure} 136 137 You should see something similar to 138 Figure~\ref{fig:Rwind2}. Surprise---the estimates coincide! To get 139 out, just close the \app{R} window or type \verb|q()| at the \app{R} 140 prompt. 141 142 \subsection{Time series data} 143 \label{sec:R-ols-arma} 144 145 We now turn to an example which uses time series data: we 146 will compare gretl's and \app{R}'s estimates of Box and Jenkins' 147 immortal airline'' model. The data are contained in the \texttt{bjg} 148 sample dataset. The following gretl code 149 \begin{code} 150 open bjg 151 arima 0 1 1 ; 0 1 1 ; lg --nc 152 \end{code} 153 produces the estimates shown in Table \ref{tab:airline-gretl}. 154 155 \begin{table}[htbp] 156 \caption{Airline model from Box and Jenkins (1976) -- selected 157 portion of gretl's estimates} 158 \label{tab:airline-gretl} 159 \begin{center} 160 161 \begin{tabular*}{0.75\textwidth}{@{\extracolsep{\fill}} 162 l% col 1: varname 163 D{.}{.}{-1}% col 2: coeff 164 D{.}{.}{-1}% col 3: sderr 165 D{.}{.}{-1}% col 4: t-stat 166 D{.}{.}{4}}% col 5: p-value (or slope) 167 Variable & 168 \multicolumn{1}{c}{Coefficient} & 169 \multicolumn{1}{c}{Std.\ Error} & 170 \multicolumn{1}{c}{t-statistic} & 171 \multicolumn{1}{c}{p-value} \\[1ex] 172 \theta_{1} & 173 -0.401824 & 174 0.0896421 & 175 -4.4825 & 176 0.0000 \\ 177 \Theta_{1} & 178 -0.556936 & 179 0.0731044 & 180 -7.6184 & 181 0.0000 \\ 182 \end{tabular*} 183 184 \begin{tabular}{lD{.}{.}{-1}} 185 Variance of innovations & 0.00134810 \\ 186 Log-likelihood & 244.696 \\ 187 Akaike information criterion & -483.39 188 \end{tabular} 189 \end{center} 190 \end{table} 191 192 If we now open an \app{R} session as described in the previous 193 subsection, the data-passing mechanism is slightly different. Since 194 our data were defined in gretl as time series, we use an \app{R} 195 \emph{time-series} object (\emph{ts} for short) for the transfer. In 196 this way we can retain in \app{R} useful information such as the 197 periodicity of the data and the sample limits. The downside is that 198 the names of individual series, as defined in gretl, are not 199 valid identifiers. In order to extract the variable \texttt{lg}, one 200 needs to use the syntax \verb|lg <- gretldata[, "lg"]|. 201 202 ARIMA estimation can be carried out by issuing the following two 203 \app{R} commands: 204 \begin{code} 205 lg <- gretldata[, "lg"] 206 arima(lg, c(0,1,1), seasonal=c(0,1,1)) 207 \end{code} 208 209 which yield 210 211 \begin{code} 212 Coefficients: 213 ma1 sma1 214 -0.4018 -0.5569 215 s.e. 0.0896 0.0731 216 217 sigma^2 estimated as 0.001348: log likelihood = 244.7, aic = -483.4 218 \end{code} 219 220 Happily, the estimates again coincide. 221 222 \section{Running an \app{R} script} 223 \label{sec:R-scripts} 224 225 Opening an \app{R} window and keying in commands is a convenient 226 method when the job is small. In some cases, however, it would be 227 preferable to have \app{R} execute a script prepared in advance. 228 One way to do this is via the \texttt{source()} command in 229 \app{R}. Alternatively, gretl offers the facility to edit an 230 \app{R} script and run it, having the current dataset pre-loaded 231 automatically. This feature can be accessed via the File, Script 232 Files'' menu entry. By selecting User file'', one can load a 233 pre-existing \app{R} script; if you want to create a new script 234 instead, select the New script, R script'' menu entry. 235 236 \begin{figure}[htbp] 237 \centering 238 \includegraphics[scale=0.7]{figures/R-edit1} 239 \caption{Editing window for \app{R} scripts} 240 \label{fig:R-edit1} 241 \end{figure} 242 In either case, you are presented with a window very similar to 243 the editor window used for ordinary gretl scripts, as in 244 Figure~\ref{fig:R-edit1}. 245 246 There are two main differences. First, you get syntax highlighting for 247 \app{R}'s syntax instead of gretl's. Second, clicking on the 248 Execute button (the gears icon), launches an instance of \app{R} in 249 which your commands are executed. Before \app{R} is actually run, you 250 are asked if you want to run \app{R} interactively or not (see 251 Figure~\ref{fig:R-exec-mode}). 252 253 \begin{figure}[htbp] 254 \centering 255 \includegraphics[scale=0.7]{figures/R-exec-mode} 256 \caption{Editing window for \app{R} scripts} 257 \label{fig:R-exec-mode} 258 \end{figure} 259 260 An interactive run opens an \app{R} instance similar to the one seen 261 in the previous section: your data will be pre-loaded (if the 262 pre-load data'' box is checked) and your commands will 263 be executed. Once this is done, you will find yourself at the \app{R} 264 prompt, where you can enter more commands. 265 266 A non-interactive run, on the other hand, will execute your script, 267 collect the output from \app{R} and present it to you in an output 268 window; \app{R} will be run in the background. If, for example, the 269 script in Figure~\ref{fig:R-edit1} is run non-interactively, a window 270 similar to Figure~\ref{fig:R-output1} will appear. 271 272 \begin{figure}[htbp] 273 \centering 274 \includegraphics[scale=0.7]{figures/R-output1} 275 \caption{Output from a non-interactive \app{R} run} 276 \label{fig:R-output1} 277 \end{figure} 278 279 \section{Taking stuff back and forth} 280 \label{sec:R-passing-data} 281 282 As regards the passing of data between the two programs, so far we 283 have only considered passing series from gretl to \app{R}. In 284 order to achieve a satisfactory degree of interoperability, more is 285 needed. In the following sub-sections we see how matrices can be 286 exchanged, and how data can be passed from \app{R} back to 287 gretl. 288 289 \subsection{Passing matrices from gretl to \app{R}} 290 291 For passing matrices from gretl to \app{R}, you can use the 292 \texttt{mwrite} matrix function described in section 293 \ref{sec:matrix-csv}. For example, the following gretl code 294 fragment generates the matrix 295 \[ 296 A = \left[ 297 \begin{array}{ccc} 298 3 & 7 & 11 \\ 299 4 & 8 & 12 \\ 300 5 & 9 & 13 \\ 301 6 & 10 & 14 302 \end{array} 303 \right] 304$
305 and stores it into the file \texttt{mymatfile.mat} in the user's
306 dotdir''(see section~\ref{sec:named-strings}). Note that writing to
307 this special directory, which is sure to exist and be writable by the
308 user, is mandated by the non-zero value for the third, optional
309 argument to \texttt{mwrite}.
310 \begin{code}
311   matrix A = mshape(seq(3,14),4,3)
312   err = mwrite(A, "mymatfile.mat", 1)
313 \end{code}
314 The recommended \app{R} code to import such a matrix is
315 \begin{code}
317 \end{code}
318
319 The function \texttt{gretl.loadmat}, which is predefined when \app{R}
320 is called from gretl, retrieves the matrix from dotdir.  (The
321 \texttt{.mat}'' extension for gretl matrix files is not compulsory;
322 you can name these files as you wish.)
323
324 It's also possible to take more control over the details of the
325 transer if you wish. You have the built-in string variable
326 \verb|$dotdir| in gretl, while in \app{R} you have the same variable 327 under the name \texttt{gretl.dotdir}. To use a location other than 328 \verb|$dotdir| you may (a) omit the third argument to \texttt{mwrite}
329 and supply a full path to the matrix file, and (b) use a more generic
330 approach to reading the file in \app{R}. Here's an example:
331
332 Gretl side:
333 \begin{code}
334   mwrite(A, "/path/to/mymatfile.mat")
335 \end{code}
336 \app{R} side:
337 \begin{code}
339 \end{code}
340
341 \subsection{Passing data from \app{R} to gretl}
342 \label{sec:Rpassing-data}
343
344 For passing data in the opposite direction, gretl defines a special
345 function that can be used in the \app{R} environment. An \app{R}
346 object will be written as a temporary file in
347 \verb|$dotdir|, from where it can be easily retrieved from within 348 gretl. 349 350 The name of this function is \texttt{gretl.export()}; it takes one 351 required argument, the object to be exported. At present, the objects 352 that can be exported with this method are matrices, data frames and 353 time-series objects. The function creates a text file, by default with 354 the same name as the exported object (plus an appropriate suffix), 355 in gretl's temporary directory. Data frames and time-series objects 356 are stored as CSV files, and can be retrieved by using gretl's 357 \texttt{append} command. Matrices are stored in a special text format 358 that is understood by gretl (see section~\ref{sec:matrix-csv}); the 359 file suffix is in this case \texttt{.mat}, and to read the matrix in 360 gretl you must use the \texttt{mread()} function. 361 362 This function also has an optional second argument, namely a string 363 which specifies a basename for the export file, in case you want to 364 use a name other than that attached to the object within \app{R}. As 365 in the default case an appropriate suffix, \texttt{.csv} or 366 \texttt{.mat}, will be added to the basename. 367 368 As an example, we take the airline data and use them to estimate a 369 structural time series model \a la \cite{harvey89}.\footnote{The 370 function package \package{StucTiSM} is available to handle this 371 class of models natively in \app{gretl}.} The model we will use is 372 the \emph{Basic Structural Model} (BSM), in which a time series is 373 decomposed into three terms: 374 $375 y_t = \mu_t + \gamma_t + \varepsilon_t 376$ 377 where$\mu_t$is a trend component,$\gamma_t$is a seasonal component 378 and$\varepsilon_t$is a noise term. In turn, the following is assumed 379 to hold: 380 \begin{eqnarray*} 381 \Delta \mu_t & = & \beta_{t-1} + \eta_t \\ 382 \Delta \beta_t & = & \zeta_t \\ 383 \Delta_s \gamma_t & = & \Delta \omega_t 384 \end{eqnarray*} 385 where$\Delta_s$is the seasonal differencing operator,$(1-L^s)$, and 386$\eta_t$,$\zeta_t$and$\omega_t$are mutually uncorrelated white 387 noise processes. The object of the analysis is to estimate the 388 variances of the noise components (which may be zero) and to recover 389 estimates of the latent processes$\mu_t$(the level''),$\beta_t$390 (the slope'') and$\gamma_t$. 391 392 We will use \app{R}'s \texttt{StructTS} command and import the results 393 back into gretl. Once the \texttt{bjg} dataset is loaded in gretl, we 394 pass the data to \app{R} and execute the following script: 395 \begin{code} 396 # extract the log series 397 y <- gretldata[, "lg"] 398 # estimate the model 399 strmod <- StructTS(y) 400 # save the fitted components (smoothed) 401 compon <- as.ts(tsSmooth(strmod)) 402 # save the estimated variances 403 vars <- as.matrix(strmod$coef)
404
405 # export into gretl's temp dir
406 gretl.export(compon)
407 gretl.export(vars)
408 \end{code}
409 %$410 411 Running this script via gretl produces minimal output: 412 \begin{code} 413 current data loaded as ts object "gretldata" 414 wrote /home/cottrell/.gretl/compon.csv 415 wrote /home/cottrell/.gretl/vars.mat 416 \end{code} 417 However, we are now able to pull the results back into gretl by 418 executing the following commands, either from the console or by 419 creating a small script: 420 \begin{code} 421 string fname = sprintf("%s/compon.csv",$dotdir)
422 append @fname
424 \end{code}
425 The first command reads the estimated time-series components from a
426 CSV file, which is the format that the passing mechanism employs for
427 series. The matrix \texttt{vars} is read from the file
428 \texttt{vars.mat}.
429
430 \begin{figure}[htbp]
431   \centering
432   \includegraphics{figures/BSM-output}
433   \caption{Estimated components from BSM}
434   \label{fig:BSM-output}
435 \end{figure}
436
437 After the above commands have been executed, three new series will
438 have appeared in the gretl workspace, namely the estimates of
439 the three components; by plotting them together with the original
440 data, you should get a graph similar to
441 Figure~\ref{fig:BSM-output}. The estimates of the variances can be
442 seen by printing the \texttt{vars} matrix, as in
443
444 \begin{code}
445 ? print vars
446 vars (4 x 1)
447
448   0.00077185
449       0.0000
450    0.0013969
451       0.0000
452 \end{code}
453
454 That is,
455 \begin{equation*}
459   \hat{\sigma}^2_{\varepsilon} = 0
460 \end{equation*}
461 Notice that, since $\hat{\sigma}^2_{\zeta} = 0$, the estimate for
462 $\beta_t$ is constant and the level component is simply a random walk
463 with a drift.
464
465 \section{Interacting with \app{R} from the command line}
466 \label{sec:foreign-command}
467
468 Up to this point we have spoken only of interaction with \app{R} via
469 the GUI program. In order to do the same from the command line
470 interface, gretl provides the \texttt{foreign} command. This
471 enables you to embed non-native commands within a gretl
472 script.
473
474 A foreign'' block takes the form
475 \begin{code}
476 foreign language=R [--send-data[=list]] [--quiet]
477     ... R commands ...
478 end foreign
479 \end{code}
480 and achieves the same effect as submitting the enclosed \app{R}
481 commands via the GUI in the non-interactive mode (see
482 section~\ref{sec:R-scripts} above). The \option{send-data} option
483 arranges for auto-loading of the data present in the gretl session, or
484 a subset thereof specified via a named list.  The \option{quiet}
485 option prevents the output from \app{R} from being echoed in the gretl
486 output.
487
488 \begin{script}[htbp]
489   \caption{Estimation of the Basic Structural Model -- simple}
490 \begin{scode}
491 open bjg.gdt
492
493 foreign language=R --send-data
494     y <- gretldata[, "lg"]
495     strmod <- StructTS(y)
496     compon <- as.ts(tsSmooth(strmod))
497     vars <- as.matrix(strmod$coef) 498 gretl.export(compon) 499 gretl.export(vars) 500 end foreign 501 502 append @dotdir/compon.csv 503 rename level lg_level 504 rename slope lg_slope 505 rename sea lg_seas 506 507 vars = mread("vars.mat", 1) 508 \end{scode} 509 \label{RStructTS-simple} 510 \end{script} 511 %$
512
513 Using this method, replicating the example in the previous subsection
514 is rather easy: basically, all it takes is encapsulating the content
515 of the \app{R} script in a \texttt{foreign}\ldots\texttt{end foreign}
516 block; see example \ref{RStructTS-simple}.
517
518 \begin{script}[htbp]
519   \caption{Estimation of the Basic Structural Model -- via a function}
520 \begin{scode}
521 function list RStructTS(series myseries)
522     smpl ok(myseries) --restrict
523     sx = argname(myseries)
524
525     foreign language=R --send-data --quiet
526         @sx <- gretldata[, "myseries"]
527         strmod <- StructTS(@sx)
528         compon <- as.ts(tsSmooth(strmod))
529         gretl.export(compon)
530     end foreign
531
532     append @dotdir/compon.csv
533     rename level @sx_level
534     rename slope @sx_slope
535     rename sea @sx_seas
536
537     list ret = @sx_level @sx_slope @sx_seas
538     return ret
539 end function
540
541 # ------------ main -------------------------
542
543 open bjg.gdt
544 list X = RStructTS(lg)
545 \end{scode}
546
547 \label{RStructTS-func}
548 \end{script}
549
550 The above syntax, despite being already quite useful by itself, shows
551 its full power when it is used in conjunction with user-written
552 functions.  Listing~\ref{RStructTS-func} shows how to define a
553 gretl function that calls \app{R} internally.
554
555 \section{Performance issues with \app{R}}
556 \label{sec:R-performance}
557
558 \app{R} is a large and complex program, which takes an appreciable
559 time to initialize itself.  In interactive use this not a significant
560 problem, but if you have a gretl script that calls \app{R} repeatedly
561 the cumulated start-up costs can become bothersome.  To get around
562 this, gretl calls the \app{R} shared library by preference; in this
563 case the start-up cost is borne only once, on the first invocation of
564 \app{R} code from within gretl.
565
566 Support for the \app{R} shared library is built into the gretl
567 packages for MS Windows and OS X---but the advantage is realized
568 only if the library is in fact available at run time.  If you are
569 building gretl yourself on Linux and wish to make use of the
570 \app{R} library, you should ensure (a) that \app{R} has been built
571 with the shared library enabled (specify \verb|--enable-R-shlib| when
572 configuring your build of \app{R}), and (b) that the \verb|pkg-config|
573 program is able to detect your \app{R} installation.  We do not link
574 to the \app{R} library at build time, rather we open it dynamically on
575 demand. The gretl GUI has an item under the
576 \textsf{Tools/Preferences} menu which enables you to select the
577 path to the library, if it is not detected automatically.
578
579 If you have the \app{R} shared library installed but want to force
580 gretl to call the \app{R} executable instead, you can do
581 \begin{code}
582 set R_lib off
583 \end{code}
584
585 \section{Further use of the \app{R} library}
586 \label{sec:R-functions}
587
588 Besides improving performance, as noted above, use of the \app{R}
589 shared library makes possible a further refinement.  That is, you can
590 define functions in \app{R}, within a \texttt{foreign} block, then
591 call those functions later in your script much as if they were
592 gretl functions.  This is illustrated below.
593 %
594 \begin{code}
595 set R_functions on
596 foreign language=R
597   plus_one <- function(q) {
598      z = q+1
599      invisible(z)
600   }
601 end foreign
602
603 scalar b=R.plus_one(2)
604 \end{code}
605 %
606 The \app{R} function \verb|plus_one| is obviously trivial in itself,
607 but the example shows a couple of points.  First, for this mechanism
608 to work you need to enable \verb|R_functions| via the \texttt{set}
609 command.  Second, to avoid collision with the gretl function
610 namespace, calls to functions defined in this way must be prefixed
611 with \texttt{R.}'', as in \verb|R.plus_one|. But note, this
612 mechanism will not work if you have defined a gretl bundle named
613 \texttt{R}: in that case identifiers beginning with \texttt{R.{}}''
614 will be understood as referring to members of the bundle in question.
615
616 Built-in \app{R} functions may also be called in this way, once
617 \verb|R_functions| is set \texttt{on}.  For example one can invoke
618 \app{R}'s \texttt{choose} function, which computes binomial
619 coefficients:
620 %
621 \begin{code}
622 set R_functions on
623 scalar b = R.choose(10,4)
624 \end{code}
625 %
626 Note, however, that the possibilities for use of built-in \app{R}
627 functions are limited; only functions whose arguments and return
628 values are sufficiently generic (basically scalars, matrices or
629 strings) will work.
630
631 %%% Local Variables:
632 %%% mode: latex
633 %%% TeX-master: "gretl-guide"
634 %%% End:
635
`