Portfolio Management with R

This book describes how to use the PMwR package. PMwR provides a small set of reliable, efficient and convenient tools that help in processing and analysing trade and portfolio data.

PMwR grew out of various pieces of software that I have written since 2008, first at the University of Geneva, later during my work at financial firms. The package has been and still is under active development, but in its current form is has become fairly stable, notably functions such as btest, returns or position. (Mainly because I have written so much code that relies on those functions that I cannot afford to change them.) Nevertheless, several parts still need some grooming, so I reserve the right to change them. Any such changes will be announced in the NEWS file (including instructions how to adapt exisiting code), and changes will be reflected in this manual. I am grateful for comments and suggestions. The book itself is written in Org mode.

The complete tangled code is available from the book's website.

The PMwR package is on CRAN; to install, simply type

install.packages("PMwR")

in an R session. The development version of the package is available from http://enricoschumann.net/R/packages/PMwR/. To get the devel version from within R, use the following command:

install.packages("PMwR",
                 repos = c('http://enricoschumann.net/R',
                           getOption('repos')))

The package depends on several other packages, which are automatically obtained from the same repository and from CRAN. The package source code is also pushed to public repositories at

https://git.sr.ht/~enricoschumann/PMwR

and

https://gitlab.com/enricoschumann/PMwR

and

https://github.com/enricoschumann/PMwR.

Recent versions of the package (since 0.3-4, released in 2016) are pure R code and can be built without any prerequisites except an R installation; older versions contained C code, so you needed to have the necessary tool chain installed (typically by installing Rtools).

I have benefited and learnt a lot from using and studying various pieces of free software.

A great inspiration were and are the principles behind the tools that make working on a Unix-type operating system so productive (and so much fun). In particular, that a programme should do one thing, and that it should do this one thing thoroughly.

Another inspiration – that small is beautiful – came from the Lua language, which I used much in the past. (Lua even removed functions.)

Closer to finance, there is the Ledger project.

small

The aim of PMwR is to provide a small set of tools. This comes at the price: interfaces may be more complicated, with arguments being overloaded. But with few functions, it is easier to remember a function name or to find it in the first place.

flexible and general

PMwR aims to be open to different types of instruments, different timestamps, etc. In this respect, I learnt much from using the zoo package. It gives a good example (I think) how a package should work: blend well with standard data structures, be idiomatic, flexible with regards to timestamps, …

functional

Start with a quote from K&R, chapter 1: "With properly designed functions, it is possible to ignore how a job is done; knowing what is done is sufficient". PMwR emphasizes computations: transforming simple and transparent data structures into other data structures. In R, a computation means calling a function. There are many good reasons for using functions.

• clearer code; easier to reuse; easier to maintain
• provide a clear view of what is needed for a specific computation; thus, they help with parallel/distributed computing
• it is easier to test functionality
• input data is not changed, which makes code more reliable
• clean workspace after function call has ended

There are more advantages, actually; such as the application of techniques such as memoisation.

Computations provided by PMwR do not – for developers: should not – rely on global options/settings. The exception are functions that are used interactively, which essentially means print methods. (In scripts or methods, cat is preferred.)

matching by name

Whenever possible and intuitive, data should be matched by name, not by position. (In other languages, we would use a hashtable instead of a vector.) This is most natural with vectors that store scalar information about instruments, such as prices or multipliers. In such cases, data input is preferred in the form of named vectors.

vectorization

Functions should do vectorization when it is beneficial in terms of speed or clarity of code. Likewise, functions should work on matrices directly (typically columnwise) when it simplifies or speeds up things. Otherwise, applying the function (i.e. looping) should be left to the user.

An example may clarify this: drawdown (actually in the NMOF package) is internally computed through cumsum, so it will be fast for a single vector. But for a matrix of time series, it would need a loop, which will be left to the user. On the other hand, returns can be computed efficiently for a matrix, so function returns directly handles matrices of prices.

Several other packages originated from PMwR; initially, much of their code had been part of PMwR. Most of those packages are on CRAN as well.

In this manual, all R output will be presented in English. In case you run R in a different locale, but want to receive messages in English, type this:

Sys.setenv(LANGUAGE = "en")

R code is shown in a typewriter font like this.

1+1

The results of a computation are shown as follows:

[1] 2

The ultimate basis of many financial computations are lists of transactions. PMwR provides an S3 class journal for handling such lists. A journal is a list of atomic vectors, to which a class attribute is attached. (Thus, a journal is similar to a dataframe.¹) Such a list is created through the function journal. Methods should not rely on this list being sorted in any particular way: components of a journal should always be retrieved by name, never by position. (In this respect a journal differs from a dataframe, for which we can meaningfully refer to the n-th column.)

A journal's components, such as amount or timestamp, are called fields in this manual.

The simplicity of the class is intended, because it is meant for interactive analyses. The user may – and is expected to – dissect the information in a journal at will; such dissections include removing the class attribute.

What is actually stored in a journal is up to the user. A number of fields are, however, required for certain operations and so it is recommended that they be present:

amount

The notional amount that is transacted. amount is, in a way, the most important property of a journal. When functions compute statistics from the journal (the number of transactions, say), they will often look at amount.

timestamp

When did the transaction take place? A numeric or character vector; should be sortable.

price

Well, price. Well, there are many types of prices. The price specified in a journal may be used to compute profit/loss; therefore, the difference between prices should be proportional to profit/loss for the transactions. Unfortunately, there are many instruments that are not quoted in transaction prices: Options may be quoted in implied volatility, bonds in yield, etc. For such instruments, the net-asset value (NAV) should be used.

instrument

Description of the financial instrument; typically an identifier, a.k.a. ticker or symbol. That is, a string or perhaps a number;² but not a more-complex object (recall that journals are lists of atomic vectors).

id

A transaction identifier, possibly but not necessarily unique.

account

Description of the account (a string or perhaps a number).

...

other fields. They must be named, as in fees = c(1, 2, 1).

All fields except amount can be missing. Such missing values will be `added back' as NA with the exception of id and account, which will be NULL. To be clear: amount could be a vector of only NA values, but it cannot be left out when the journal is created. (This will become clearer with the examples below.)

A journal may have no transactions at all in it. In such a case all fields have length zero,

e.g. amount would be numeric(0) and so on. Such empty journals may be created by saying journal() or by coercing a zero-row dataframe to a journal, via a call to as.journal.

Transactions in a journal may be organised in hierarchies, such as

account => subaccount => subsubaccont => ... => instrument

This is useful and necessary when you have traded an instrument for different accounts, say, or as part of different strategies. Such a hierarchy may be completely captured in the instrument field, by concatenating account hierarchy and instrument using a separator pattern such as ::.³ The result would be 'namespaced' instruments such as Pension::Equities::AMZN. Alternatively, part of the hierarchy may be stored in the account field.

The function journal creates journal objects. See ?journal for details about the function and about methods for journal objects. At its very minimum, a journal must contain amounts of something.

J <- journal(amount = c(1, 2, -2, 3))
J

   amount
1       1
2       2
3      -2
4       3

4 transactions

Actually, that is not true. Sometimes it is useful to create an empty journal, one with no entries at all. You can do so by saying journal(), without any arguments.

journal()

no transactions

To see the current balance, which is nothing more than the sum over all amounts, you can use position.

position(J)

4

Suppose you wanted to note how many bottles of milk and wine you have stored in your basement. Whenever you add to your stock, you have a positive amount; whenever you retrieve bottles, you have a negative amount. Then, by keeping track of transactions, you do not have to take stock (apart, perhaps, from occasional checking that you did not miss a transaction), as long as you keep track of what you put into your cellar and what you take out.

There may be some analyses you can do on flows alone: you may check your drinking habits for patterns, such as slow accumulation of wine, followed by rapid consumption; or the other way around. But typically, you will want to analyse your transactions later, and then the more information you record about them – when, what, why, at what price, etc. –, the better. Journals allow you to store such information. To show how they are used, let us switch to a financial example.

Suppose you have transacted the following trades.

|  timestamp | account | instrument | amount |   price | note        |
|------------+---------+------------+--------+---------+-------------|
| 2017-08-01 | Pension | AMZN       |     10 | 1001.00 |             |
| 2017-08-01 | Pension | MSFT       |    220 |   73.10 |             |
| 2017-07-14 | Trading | AMZN       |     10 | 1001.50 |             |
| 2017-07-31 | Trading | AMZN       |     -5 | 1014.00 | take profit |
| 2017-08-15 | Trading | AMZN       |     10 |  985.50 |             |
| 2017-10-05 | Pension | MSFT       |     70 |   74.40 |             |

This table is formatted in Org syntax. The function \iRfunction{org\_journal} converts such a table into a journal.

org_journal <- function(file, text, timestamp.as.Date = TRUE) {
  ans <- orgutils::readOrg(text = text)
  if (timestamp.as.Date && "timestamp" %in% colnames(ans))
    ans$timestamp <- as.Date(ans$timestamp)
  ans <- as.journal(ans)
  ans
}

We read the table into a journal.

J <- org_journal(text = "
  |  timestamp | account | instrument | amount |   price | note        |
  |------------+---------+------------+--------+---------+-------------|
  | 2017-08-01 | Pension | AMZN       |     10 | 1001.00 |             |
  | 2017-08-01 | Pension | MSFT       |    220 |   73.10 |             |
  | 2017-07-14 | Trading | AMZN       |     10 | 1001.50 |             |
  | 2017-07-31 | Trading | AMZN       |     -5 | 1014.00 | take profit |
  | 2017-08-15 | Trading | AMZN       |     10 |  985.50 |             |
  | 2017-10-05 | Pension | MSFT       |     70 |   74.40 |             |
    ")
J

A print method defines how a journal is displayed. See ?print.journal for details. In general, you can always get help for methods for generic functions by saying ?<generic_function>.journal, e.g. ?print.journal or ?as.data.frame.journal.

print(J, max.print = 2, exclude = c("account", "note"))

   instrument   timestamp  amount   price
1        AMZN  2017-08-01      10  1001.0
2        MSFT  2017-08-01     220    73.1
[ .... ]

6 transactions

A str method shows the fields in the journal.

str(J)

'journal':	 6 transactions
 $ instrument: chr [1:6] "AMZN" "MSFT" "AMZN" "AMZN" ...
 $ account   : chr [1:6] "Pension" "Pension" "Trading" ...
 $ timestamp : Date[1:6], format: "2017-08-01" ...
 $ amount    : int [1:6] 10 220 10 -5 10 70
 $ price     : num [1:6] 1001 73.1 1001.5 1014 985.5 ...
 $ note      : chr [1:6] "" "" "" "take profit" ...

You may notice that the output is similar to that of a data.frame or list. That is because J is a list of atomic vectors, with a class attribute. Essentially, it is little more than a list of the columns of the above table.

But note that journal would silently have added required fields such price, initialised as NA.

str(journal(amount = 1))

'journal':     1 transaction
 $ instrument: chr NA
 $ timestamp : num NA
 $ amount    : num 1
 $ price     : num NA

In the example, the timestamps are of class Date. But essentially, any vector of mode character or numeric can be used, for instance POSIXct, or other classes. Here is an example that uses the \iRpackage{nanotime} package \citep{eddelbuettel2017}.

library("nanotime")
journal(amount = 1:3, timestamp = nanotime(Sys.time()) + 1:3)

                             timestamp  amount
1  2020-08-09T08:27:19.812951001+00:00       1
2  2020-08-09T08:27:19.812951002+00:00       2
3  2020-08-09T08:27:19.812951003+00:00       3

3 transactions

Journals can be combined with c.

J2 <- J
J2$remark <- rep("new", length(J))
c(J, J2)

    instrument   timestamp  amount   price  account         note  remark
1         AMZN  2017-08-01      10  1001.0  Pension                 <NA>
2         MSFT  2017-08-01     220    73.1  Pension                 <NA>
3         AMZN  2017-07-14      10  1001.5  Trading                 <NA>
4         AMZN  2017-07-31      -5  1014.0  Trading  take profit    <NA>
5         AMZN  2017-08-15      10   985.5  Trading                 <NA>
6         MSFT  2017-10-05      70    74.4  Pension                 <NA>
7         AMZN  2017-08-01      10  1001.0  Pension                  new
8         MSFT  2017-08-01     220    73.1  Pension                  new
9         AMZN  2017-07-14      10  1001.5  Trading                  new
10        AMZN  2017-07-31      -5  1014.0  Trading  take profit     new
11        AMZN  2017-08-15      10   985.5  Trading                  new
12        MSFT  2017-10-05      70    74.4  Pension                  new

12 transactions

The new combined journal will not be sorted by date. In general, a journal need not be sorted in any particular way. There is a sort method for journals, whose default is to sort by timestamp.

We can also sort by other fields, for instance by amount.

sort(c(J, J2), by = c("amount", "price"), decreasing = FALSE)

    instrument   timestamp  amount   price  account         note  remark
1         AMZN  2017-07-31      -5  1014.0  Trading  take profit    <NA>
2         AMZN  2017-07-31      -5  1014.0  Trading  take profit     new
3         AMZN  2017-08-15      10   985.5  Trading                 <NA>
4         AMZN  2017-08-15      10   985.5  Trading                  new
5         AMZN  2017-08-01      10  1001.0  Pension                 <NA>
6         AMZN  2017-08-01      10  1001.0  Pension                  new
7         AMZN  2017-07-14      10  1001.5  Trading                 <NA>
8         AMZN  2017-07-14      10  1001.5  Trading                  new
9         MSFT  2017-10-05      70    74.4  Pension                 <NA>
10        MSFT  2017-10-05      70    74.4  Pension                  new
11        MSFT  2017-08-01     220    73.1  Pension                 <NA>
12        MSFT  2017-08-01     220    73.1  Pension                  new

12 transactions

You can query the number of transactions in a journal with length.

In an interactive session, you can use subset to select transactions.

subset(J, amount > 10)

   instrument   timestamp  amount  price  account  comment
1        MSFT  2017-08-01     220   73.1  Pension
2        MSFT  2017-10-05      70   74.4  Pension

2 transactions

With subset, you need not quote the expression that selects trades and you can directly access a journal's fields. Because of the way subset evaluates its arguments, it should not be used within functions. (See the Examples section in ?journal for what can happen then.)

More generally, to extract or change a field, use its name, either through the $ operator or double brackets [[...]].⁴

J$amount

[1]  10 220  10  -5  10  70

You can also replace specific fields.

J$comment[1] <-  "some note"
J

   instrument   timestamp  amount   price  account      comment
1        AMZN  2017-08-01      10  1001.0  Pension       a note
2        MSFT  2017-08-01     220    73.1  Pension
3        AMZN  2017-07-14      10  1001.5  Trading
4        AMZN  2017-07-31      -5  1014.0  Trading  take profit
5        AMZN  2017-08-15      10   985.5  Trading
6        MSFT  2017-10-05      70    74.4  Pension

6 transactions

The `[` method works with integers or logicals, returning the respective transactions.

J[2:3]

   instrument   timestamp  amount   price  account  comment
1        MSFT  2017-08-01     220    73.1  Pension
2        AMZN  2017-07-14      10  1001.5  Trading

2 transactions

J[J$amount < 0]

   instrument   timestamp  amount  price  account      comment
1        AMZN  2017-07-31      -5   1014  Trading  take profit

1 transaction

You can also pass a string, which is then interpreted as a regular expression that is matched against all character fields in the journal.

J["Pension"]

   instrument   timestamp  amount   price  account  comment
1        AMZN  2017-08-01      10  1001.0  Pension   a note
2        MSFT  2017-08-01     220    73.1  Pension
3        MSFT  2017-10-05      70    74.4  Pension

3 transactions

You can also specify the fields to match the string against.

J["Pension", match.against = "instrument"]

no transactions

By default, case is ignored, but you can set ignore.case to FALSE.

J["pension", ignore.case = FALSE]

no transactions

Finally, you can invert the selection with invert.

J["Pension", invert = TRUE]

   instrument   timestamp  amount   price  account      comment
1        AMZN  2017-07-14      10  1001.5  Trading
2        AMZN  2017-07-31      -5  1014.0  Trading  take profit
3        AMZN  2017-08-15      10   985.5  Trading

3 transactions

Often the data provided by journals needs to be processed in some way. A straightforward strategy is to call as.data.frame on the journal and then to use one of the many functions and methods that can be used for dataframes, such as aggregate or apply.

Even without coercion to a dataframe: A journal is a list of atomic vectors and hence already very similar to a dataframe. As a consequence, many computations can also be done directly on the journal, in particular with tapply. An example: you have a journal trades and want to compute monthly turnover (two-way). If there is only one instrument or all instruments may be added without harm (typically when they are denominated in the same currency), you can use this expression:

tapply(trades,
       INDEX = format(trades$timestamp, "%Y-%m"),
       FUN = function(x) sum(abs(x$amount)))

To break it down by instrument, just add instrument as a second grouping variable to the INDEX argument.

tapply(trades,
       INDEX = list(format(trades$timestamp, "%Y-%m"),
                    trades$instrument),
       FUN = function(x) sum(abs(x$amount)))

A special case is when a journal is to be processed into a new journal. For this, PMwR defines an aggregate method for journals:

aggregate.journal(x, by, FUN, ...)

The method splits the journal according to the grouping argument by, which can be a list (as in the default method) or an atomic vector. The argument FUN can either be a function or list. If a function, it should receive a journal and also evaluate to a journal. (Note that this is different from R's aggregate.data.frame, which calls FUN on all columns, but in turn cannot address specific columns of the data.frame.) If FUN is a list, its elements should be named functions. The names should match fields in the journal.

An example: we have a journal covering two trading days and wish to create a summary journal, which aggregates buys and sells for every day.

J <- org_journal(text = "
  | instrument | timestamp      | amount | price |
  |------------+----------------+--------+-------|
  | A          | 2013-09-02 Mon |     -3 |   102 |
  | B          | 2013-09-02 Mon |     -3 |   104 |
  | B          | 2013-09-02 Mon |      3 |   106 |
  | B          | 2013-09-02 Mon |     -2 |   104 |
  | A          | 2013-09-03 Tue |     -1 |   110 |
  | A          | 2013-09-03 Tue |      1 |   104 |
  | A          | 2013-09-03 Tue |      5 |   108 |
  | A          | 2013-09-03 Tue |      3 |   107 |
  | B          | 2013-09-03 Tue |     -4 |   102 |
  | B          | 2013-09-03 Tue |      3 |   106 |
")

fun <- function(x) {
  journal(timestamp = as.Date(x$timestamp[1]),
          amount = sum(x$amount),
          price = sum(x$amount*x$price)/sum(x$amount),
          instrument = x$instrument[1L])
}

aggregate(J,
          by = list(J$instrument,
                    sign(J$amount),
                    as.Date(J$timestamp)),
          FUN = fun)

The results is a journal, but with at most a single buy or sell transaction per instrument per day: see the buy transaction for instrument A on September, 3.

   instrument   timestamp  amount     price
1           A  2013-09-02      -3  102.0000
2           B  2013-09-02      -5  104.0000
3           B  2013-09-02       3  106.0000
4           A  2013-09-03      -1  110.0000
5           B  2013-09-03      -4  102.0000
6           A  2013-09-03       9  107.2222
7           B  2013-09-03       3  106.0000

7 transactions

Unfortunately, in real life computing P/L is often more complicated:

• One asset-price unit may not translate into one currency unit: there may be multipliers a.k.a. contract factors; there are even instruments with variable multipliers, e.g. Australian government-bond futures.
• Asset positions may map into cashflows in non-obvious ways. The simple case is the delay in actual payment and delivery of an asset, which is often two or three days. The more problematic cases are derivatives with daily adjustments of margins. In such cases, one may need to model (i.e. keep track of) the actual account balances.
• Assets may be denominated in various currencies.
• Currencies themselves may be assets in the portfolio. Depending on how they are traded (cash, forwards, &c.), computing P/L may not be straightforward.

How – or, rather, to what degree – these complications are handled is, as always, up to the user. For a single instrument, computing P/L in units of the instrument is usually meaningful, though perhaps not always intuitive. But adding up the profits and losses of several assets will in general not work because of multipliers or different currencies. A simple and transparent way is then to manipulate the journal before P/L is computed (e.g., multiply notionals by their multipliers).

The function returns computes returns from prices. The function computes what are called simple or discrete returns:⁶ let \(P_t\) be the price at point in time \(t\), then

For computing profit/loss in currency units, see Section Computing profit and (or) loss.

Typically, we transform a whole series \(P_{t_1}, P_{t_2}, P_{t_3}, \ldots\) into returns \(R_{t_2}, R_{t_3}, \ldots\), which is a one-liner in R:

simple_returns <- function(x)
    x[-1L]/x[-length(x)] - 1

(You may argue that these are two lines: yet even a one-liner, if used repeatedly, should be written as a function.)

Let us try it. PMwR comes with two small datasets, DAX and REXP. DAX stands for Deutscher Aktienindex (German Equity Index), and REXP stands for Rentenindex (Performance). Both datasets are dataframes of one column that contains the price for the day, with the timestamps stored as rownames in format YYYY-MM-DD.

head(DAX)

               DAX
2014-01-02 9400.04
2014-01-03 9435.15
2014-01-06 9428.00
2014-01-07 9506.20
2014-01-08 9497.84
2014-01-09 9421.61

We extract the prices for the first five business days of 2014 and put them into a vector P.

P <- head(DAX[[1]], n = 5)
P

[1] 9400.04 9435.15 9428.00 9506.20 9497.84

Now we call simple_returns.

simple_returns(P)

[1]  0.00373509 -0.00075780  0.00829444 -0.00087943

In fact, using returns as provided by PMwR would have given the same result.

returns(P)

[1]  0.00373509 -0.00075780  0.00829444 -0.00087943

PMwR's returns function offers more convenience than simple_returns. For instance, it will recognise when the input argument has several columns, such as a matrix or a dataframe. In such a case, it computes returns for each column.⁷

returns(cbind(P, P))

             P         P
[1,]  0.003735  0.003735
[2,] -0.000758 -0.000758
[3,]  0.008294  0.008294
[4,] -0.000879 -0.000879

The argument pad determines how the initial observation is handled. The default, NULL, means that the first observation is dropped. It is often useful to use NA instead, since in this way the returns series keeps the same length as the original price series.

data.frame(price = P, returns = returns(P, pad = NA))

   price     returns
1 9400.0          NA
2 9435.1  0.00373509
3 9428.0 -0.00075780
4 9506.2  0.00829444
5 9497.8 -0.00087943

Setting pad to 0 can also be useful, because then it is easy to 'rebuild' the original series with cumprod. (But see Section Scaling series for a description of the function scale1, which is even more convenient.)

P[1] * cumprod(1 + returns(P, pad = 0))

[1] 9400.04 9435.15 9428.00 9506.20 9497.84

returns also has an argument lag, with default 1. This can be used to compute rolling returns, such as 30-day returns, etc.

returns is a generic function, which goes along with some overhead. If you need to compute returns on simple data structures as in the examples above and need fast computation, then you may also use .returns. The function is the actual workhorse that performs the raw return calculation.

Besides having methods for numeric vectors (which includes those with a dim attribute, i.e. matrices) and dataframes, returns also understands zoo objects. So let us create two zoo series, dax and rex.

library("zoo")
dax <- zoo(DAX[[1]], as.Date(row.names(DAX)))
rex <- zoo(REXP[[1]], as.Date(row.names(REXP)))

str(dax)

'zoo' series from 2014-01-02 to 2015-12-30
  Data: num [1:505] 9400 9435 9428 9506 9498 ...
  Index:  Date[1:505], format: "2014-01-02" "2014-01-03" ...

str(rex)

'zoo' series from 2014-01-02 to 2015-12-30
  Data: num [1:502] 441 441 442 442 442 ...
  Index:  Date[1:502], format: "2014-01-02" "2014-01-03" ...

returns(head(dax, 5), pad = NA)

2014-01-02 2014-01-03 2014-01-06 2014-01-07 2014-01-08
        NA   0.003735  -0.000758   0.008294  -0.000879

Matrices work as well. We combine both series into a two-column matrix drax.⁸

drax <- cbind(dax, rex)
returns(head(drax, 5))

                 dax       rex
2014-01-03  0.003735  0.000611
2014-01-06 -0.000758  0.001704
2014-01-07  0.008294  0.000621
2014-01-08 -0.000879 -0.000131

As you see, just as for a numeric matrix, the function computes the returns for each column.

In fact, zoo objects bring another piece of information – timestamps – that returns can use. (Since xts series inherit from zoo, they will work as well.)

When a timestamp is available, returns can compute returns for specific calendar periods. As an example, we use the daily DAX levels in 2014 and 2015 and compute monthly returns from them.

returns(dax, period = "month")

      Jan Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct Nov  Dec YTD
2014 -1.0 4.1 -1.4  0.5  3.5 -1.1 -4.3  0.7  0.0 -1.6 7.0 -1.8 4.3
2015  9.1 6.6  5.0 -4.3 -0.4 -4.1  3.3 -9.3 -5.8 12.3 4.9 -5.6 9.6

If you prefer to not use zoo or xts, you may also pass the timestamp explicitly to returns.

returns(coredata(dax), t = index(dax), period = "month")

      Jan Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct Nov  Dec YTD
2014 -1.0 4.1 -1.4  0.5  3.5 -1.1 -4.3  0.7  0.0 -1.6 7.0 -1.8 4.3
2015  9.1 6.6  5.0 -4.3 -0.4 -4.1  3.3 -9.3 -5.8 12.3 4.9 -5.6 9.6

Despite the way these monthly returns are printed: the result of the function call is a numeric vector (the return numbers), with additional information added through attributes. There is also a class attribute, which has value p_returns. The advantage of this data structure is that it is `natural' to compute with the returns, e.g. to calculate means, extremes and other quantities.

range(returns(dax, period = "month"))

[1] -0.0928  0.1232

Most useful, however, is probably the print method,

whose results you have seen above. You may also compute monthly returns for matrices, i.e. for more than one asset. But now the print method will behave differently. The function's assumption is that now it would be more convenient to print the returns aligned by date in a table.

returns(drax, period = "month")

             dax   rex
2014-01-31  -1.0   1.8
2014-02-28   4.1   0.4
2014-03-31  -1.4   0.1
2014-04-30   0.5   0.3
2014-05-30   3.5   0.9
2014-06-30  -1.1   0.4
2014-07-31  -4.3   0.4
2014-08-29   0.7   1.0
2014-09-30   0.0  -0.1
2014-10-31  -1.6   0.1
2014-11-28   7.0   0.4
2014-12-30  -1.8   1.0
2015-01-30   9.1   0.3
2015-02-27   6.6   0.1
2015-03-31   5.0   0.3
2015-04-30  -4.3  -0.5
2015-05-29  -0.4  -0.2
2015-06-30  -4.1  -0.8
2015-07-31   3.3   0.7
2015-08-31  -9.3   0.0
2015-09-30  -5.8   0.4
2015-10-30  12.3   0.4
2015-11-30   4.9   0.3
2015-12-30  -5.6  -0.6

If you rather wanted the other, one-row-per-year display, just call the function separately for each series.

lapply(list(DAX = dax, REXP = rex),
       returns, period = "month")

$DAX
      Jan Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct Nov  Dec YTD
2014 -1.0 4.1 -1.4  0.5  3.5 -1.1 -4.3  0.7  0.0 -1.6 7.0 -1.8 4.3
2015  9.1 6.6  5.0 -4.3 -0.4 -4.1  3.3 -9.3 -5.8 12.3 4.9 -5.6 9.6

$REXP
     Jan Feb Mar  Apr  May  Jun Jul Aug  Sep Oct Nov  Dec YTD
2014 1.8 0.4 0.1  0.3  0.9  0.4 0.4 1.0 -0.1 0.1 0.4  1.0 7.1
2015 0.3 0.1 0.3 -0.5 -0.2 -0.8 0.7 0.0  0.4 0.4 0.3 -0.6 0.5

See ?print.p_returns for more display options. For instance:

print(returns(dax, period = "month"),
      digits = 2, year.rows = FALSE, plus = TRUE,
      month.names = 1:12)

      2014    2015
1    -1.00   +9.06
2    +4.14   +6.61
3    -1.40   +4.95
4    +0.50   -4.28
5    +3.54   -0.35
6    -1.11   -4.11
7    -4.33   +3.33
8    +0.67   -9.28
9    +0.04   -5.84
10   -1.56  +12.32
11   +7.01   +4.90
12   -1.76   -5.62
YTD  +4.31   +9.56

There are methods toLatex and toHTML for monthly returns. In Sweave documents, you need to use results = tex and echo = false in the chunk options:

\begin{tabular}{rrrrrrrrrrrrrr}
<<results=tex,echo=false>>=
toLatex(returns(dax, period = "month"))
\end{tabular}

(There is also a vignette that gives examples for toLatex; say vignette("FinTeX", package = "PMwR") to open it.)

To get annualised returns, use period ann (or actually any string matched by the regular expression ^ann; case is ignored).

returns(dax, period = "ann")

6.9%  [02 Jan 2014 -- 30 Dec 2015]

Now let us try a shorter period.

returns(window(dax, end = as.Date("2014-1-31")),
        period = "ann")

-1.0%  [02 Jan 2014 -- 31 Jan 2014;
        less than one year, not annualised]

The function did not annualise: it refuses to do so if the time period is shorter than one year. (You may verify the return for January 2014 in the tables above.) To force annualising, add a !. The exclamation mark serves as a mnenomic that it is now imperative to annualise.

returns(window(dax, end = as.Date("2014-1-31")),
        period = "ann!")

-11.8%  [02 Jan 2014 -- 31 Jan 2014;
         less than one year, but annualised]

There are several more accepted values for period, such as year, quarter, month-to-date (mtd), year-to-date (ytd) or inception-to-date (total). The help page of returns lists all options. Note that any such setting for period requires that the timestamp can be coerced to Date; for instance, intraday time-series with POSIXct timestamps would work as well.

Sometimes we may need to compute returns for a portfolio of fixed weights, given an assumption when the portfolio is rebalanced. For instance, we may want to see how a constant allocation of 10%, 50% and 40%. to three funds would have done, assuming that a portfolio is rebalanced once a month. If more detail is necessary, then function btest can be used; see Chapter Backtesting. But the simple case can be done with returns already. Here is an example.

prices <- c(100, 102, 104, 104, 104.5,
              2, 2.2, 2.4, 2.3,   2.5,
            3.5,   3, 3.1, 3.2,   3.1)

dim(prices) <- c(5, 3)
prices

      [,1] [,2] [,3]
[1,] 100.0  2.0  3.5
[2,] 102.0  2.2  3.0
[3,] 104.0  2.4  3.1
[4,] 104.0  2.3  3.2
[5,] 104.5  2.5  3.1

Now suppose we want a constant weight vector, \([0.1, 0.5, 0.4]'\), but only rebalance at times 1 and 4. (That is, we rebalance the portfolio only with the prices at timestamps 1 and 4.)

returns(prices,
        weights = c(10, 50, 40)/100,
        rebalance.when = c(1, 4))

[1] -0.0051429  0.0637565 -0.0128240  0.0314590
attr(,"holdings")
           [,1]    [,2]    [,3]
[1,] 0.00100000 0.25000 0.11429
[2,] 0.00100000 0.25000 0.11429
[3,] 0.00100000 0.25000 0.11429
[4,] 0.00096154 0.21739 0.12500
[5,] 0.00096154 0.21739 0.12500
attr(,"contributions")
           [,1]      [,2]      [,3]
[1,] 0.00200000  0.050000 -0.057143
[2,] 0.00201034  0.050258  0.011488
[3,] 0.00000000 -0.023623  0.010799
[4,] 0.00048077  0.043478 -0.012500

The result also contains, as attributes, the imputed holdings and the single period contributions.

Argument weights does not have to be a vector. It can also be a matrix. In such a case, each row is interpreted as a portfolio. Instead of weights, we could also specify fixed positions. See ?returns for different possibilities to call returns.

Let \(w(t,i)\) be the weight of portfolio segment \(i\) at the beginning of period \(t\), and let \(r(t,i)\) be the return of segment \(i\) over period \(t\). Then the portfolio return over period \(t\), \(r_{\mbox{\scriptsize P}}(t)\) is a weighted sum of the \(N\) segment returns.

When the weights sum to unity, we may also write

or, defining \(1+r \equiv R\),

The total return contribution of segment \(i\) over time equals

In this way, a segment's return contribution in one period is reinvested in the overall portfolio in succeeding periods.

The calculation is provided in the function rc (`return contribution').

weights <- rbind(c( 0.25, 0.75),  ## the assets' weights
                 c( 0.40, 0.60),  ## during three periods
                 c( 0.25, 0.75))

R <- rbind(c( 1  ,    0),         ## the assets' returns
           c( 2.5, -1.0),         ## during these periods
           c(-2  ,  0.5))/100

rc(R, weights, segment = c("equities", "bonds"))

$period_contributions
  timestamp equities    bonds    total
1         1   0.0025  0.00000  0.00250
2         2   0.0100 -0.00600  0.00400
3         3  -0.0050  0.00375 -0.00125

$total_contributions
equities    bonds    total
 0.00749 -0.00224  0.00525

External cashflows (or transfers of positions) can be handled just like dividends. The followwing table shows the values and cashflows of a hypothetical portfolio.

A total-return series, based on column value with cashflow but excluding column cashflow, can be computed with div\_adjust.

cf <- c(100, 100, -200)
t <- c(1, 4, 5)
x <- c(100, 101, 104, 203, 4)
div_adjust(x, t, div = -cf, backward = FALSE)

[1] 100.000 101.000 104.000 103.000 103.507

More conveniently, the function unit_prices helps to compute so-called time-weighted returns of a portfolio when there are in- and outflows.

(The term time-weighted returns is actually a misnomer, as returns are not weighted at all. They are only time-weighted if time-periods are of equal length.) We repeat the previous example.

NAV <- data.frame(timestamp = 1:5,
                  NAV = x)

cf <- data.frame(timestamp = t,
                 cashflow = cf)
unit_prices(NAV, cf)

  timestamp NAV   price     units
1         1 100 100.000 1.0000000
2         2 101 101.000 1.0000000
3         3 104 104.000 1.0000000
4         4 203 103.000 1.9708738
5         5   4 103.507 0.0386446

The function returns a dataframe: to compute returns, use the price column.

At a given instant in time (in actual life, `now'), a trader needs to answer the following questions:

1. Do I want to compute a new target portfolio, yes or no? If yes, go ahead and compute the new target portfolio.
2. Given the target portfolio and the actual portfolio, do I want to rebalance (i.e. close the gap between the actual portfolio and the target portfolio)? If yes, rebalance.

If such a decision is not just hypothetical, then the answer to the second question may lead to a number of orders sent to a broker. Note that many traders do not think in terms of stock (i.e. balances) as we did here; rather, they think in terms of flow (i.e. orders). Both approaches are equivalent, but the described one makes it easier to handle missed trades and synchronise accounts.

During a backtest, we will simulate the decisions of the trader. How precisely we simulate depends on the trading strategy. The btest function is meant as a helper function to simulate these decisions. The logic for the decisions described above must be coded in the functions do.signal, signal and do.rebalance.

Implementing btest required a number of decision too: (i) what to model (i.e. how to simulate the trader), and (ii) how to code it. As an example for point (i): how precisely do we want to model the order process (e.g. use limit orders?, allow partial fills?) Example for (ii): the backbone of btest is a loop that runs through the data. Loops are slow in R when compared with compiled languages,⁹ so should we vectorise instead? Vectorisation is indeed often possible, namely if trading is not path-dependent. If we have already a list of trades, we can efficiently transform them into a profit-and-loss in R without relying on an explicit loop (see Section Computing profit and (or) loss). Yet, one advantage of looping is that the trade logic is more similar to actual trading; we may even be able to reuse some code in live trading.

Altogether, the aim for btest is to stick to the functional paradigm as much as possible. Functions receive arguments and evaluate to results; but they do not change their arguments, nor do they assign or change other variables `outside' their environment, nor do the results depend on some variable outside the function. This creates a problem, namely how to keep track of state. If we know what variables need to be persistent, we could pass them to the function and always have them returned. But we would like to be more flexible, so we can pass an environment; examples are below. To make that clear: functional programming should not be seen as a yes-or-no decision; it is a matter of degree. And more of the functional approach can help already.

All computations of btest will be based on one or several price series of length T. Internally, these prices are stored in numeric matrices.

Prices are passed as argument prices. For a single asset, this must be a matrix of prices with four columns: open, high, low and close.

For n assets, you need to pass a list of length four: prices[[1]] must be a matrix with n columns containing the open prices for the assets; prices[[2]] is a matrix with the high prices, and so on. For instance, with two assets, you need four matrices with two columns each:

 open     high     low     close
+-+-+    +-+-+    +-+-+    +-+-+
| | |    | | |    | | |    | | |
| | |    | | |    | | |    | | |
| | |    | | |    | | |    | | |
| | |    | | |    | | |    | | |
| | |    | | |    | | |    | | |
+-+-+    +-+-+    +-+-+    +-+-+

If only close prices are used, then for a single asset, use either a matrix of one column or a numeric vector. For multiple assets a list of length one must be passed, containing a matrix of close prices. For example, with 100 close prices of 5 assets, the prices should be arranged in a matrix p of size 100 times 5; and prices = list(p).

The btest function runs from b+1 to T. The variable b is the burn-in and it needs

to be a positive integer. When we take decisions that are based on past data, we will lose at least one data point. In rare cases b may be zero.

Here is an important default: at time =t=, we can use information up to time t-1. Suppose that t were 4. We may use all information up to time 3, and trade at the open in period 4:

t    time      open  high  low   close
1    HH:MM:SS                             <--\
2    HH:MM:SS                             <-- - use information
3    HH:MM:SS  _________________________  <--/
4    HH:MM:SS     X                       <- trade here
5    HH:MM:SS

We could also trade at the close:

t    time      open  high  low   close
1    HH:MM:SS                             <-- \
2    HH:MM:SS                             <-- - use information
3    HH:MM:SS  _________________________  <-- /
4    HH:MM:SS                        X    <-- trade here
5    HH:MM:SS

No, we cannot trade at the high or low. (Some people like the idea, as a robustness check, to always buy at the high, sell at the low. Robustness checks – forcing a bit of bad luck into the simulation – are a good idea, notably bad executions. High/low ranges can inform such checks, but using these ranges does not go far enough, and is more of a good story than a meaningful test.)

It is best to describe the btest function through a number of simple examples.

It does not really make a difference whether btest is called with a single or with several instruments. The pattern in signal is still to call Close() and friends to obtain data, but now these functions will return matrices with more than one column. For instance, when you have 5 assets, then Close(n = 250) would return a matrix of size 250 × 5. When signal has finished its computations, it is now expected to return a vector of positions or weights. In the example with 5  assets, it should return a vector of length 5.

There is more than one way to accomplish a certain task.

The function rebalance computes the transactions necessary for moving from one portfolio to another. The default setting is that the current portfolio is in currency units; the target portfolio in weights.

To compute the required order sizes, we also need the current prices of the assets. When current, target and price are unnamed, the assets' positions in the vectors need to match.

Suppose we have three stocks A, B and C with prices 1, 2 and 3. The main use case is a situation like this: you hold 50, 30 and 20 shares of these three stocks. However, suppose you have a target weight of 50%, 30% and 20%.

prices <- 1:3
current <- c(50, 30, 20)
target <- c(0.5, 0.3, 0.2)
rebalance(current, target, prices, match.names = FALSE)

  price current value    %       target value    %       order
1     1      50    50 29.4           85    85 50.0          35
2     2      30    60 35.3           26    52 30.6          -4
3     3      20    60 35.3           11    33 19.4          -9

Notional: 170.  Target net amount : 170.  Turnover (2-way): 70.

Or perhaps you prefer an equal weight for every asset. Note that target now is a single number.

rebalance(current, target = 1/length(current),
          price = prices, match.names = FALSE)

  price current value    %       target value    %       order
1     1      50    50 29.4           57    57 33.5           7
2     2      30    60 35.3           28    56 32.9          -2
3     3      20    60 35.3           19    57 33.5          -1

Notional: 170.  Target net amount : 170.  Turnover (2-way): 14.

Note that the target weights cannot be reached exactly because the function rounds to integers.

rebalance also supports a number of special cases. Suppose you want to go into cash and close every position.

1: rebalance(current = current, target = 0,
2:           price = prices, match.names = FALSE)

  price current value    %       target value   %       order
1     1      50    50 29.4            0     0 0.0         -50
2     2      30    60 35.3            0     0 0.0         -30
3     3      20    60 35.3            0     0 0.0         -20

Notional: 170.  Target net amount : 0.  Turnover (2-way): 170.

Suppose we have no current position and want to give equal weight to each stock. Note first that we need to specify a notional. Also, rebalance now assumes that you want to invest in every stock for which a price is supplied.

1: rebalance(current = 0, target = 1/3, notional = 100,
2:           price = prices, match.names = FALSE)

  price current value   %       target value    %       order
1     1       0     0 0.0           33    33 33.0          33
2     2       0     0 0.0           17    34 34.0          17
3     3       0     0 0.0           11    33 33.0          11

Notional: 100.  Target net amount : 100.  Turnover (2-way): 100.

More usefully, rebalance can also use the names of the vectors current, target and price. The argument match.names must be set to TRUE for this (which is the default, actually).

prices  <- c(1,1,1,1)
names(prices) <- letters[1:4]
current <- c(a = 0, b = 10)
target  <- c(a = 0, d = 0.5)
rebalance(current, target, prices)

  price current value     %       target value    %       order
b     1      10    10 100.0            0     0  0.0         -10
d     1       0     0   0.0            5     5 50.0           5

Notional: 10.  Target net amount : 5.  Turnover (2-way): 15.

To also show all instruments, set the argument drop.zero to FALSE.

print(rebalance(current, target, prices), drop.zero = FALSE)

  price current value     %       target value    %       order
a     1       0     0   0.0            0     0  0.0           0
b     1      10    10 100.0            0     0  0.0         -10
d     1       0     0   0.0            5     5 50.0           5

Notional: 10.  Target net amount : 5.  Turnover (2-way): 15.

In Section Keeping track of transactions: journals we used the function position to compute balances from transactions. The function may also be directly used to set up a position.

position(amount = 1)

1

position(amount = c(1, 2, 3), instrument = letters[1:3])

a 1
b 2
c 3

Note that with more than one instrument, these instruments need to be named. Otherwise, position would aggregate the positions (as it does for a journal). As a short-cut, you can also pass a named vector. (See http://enricoschumann.net/notes/computing-positions.html for more details.)

position(amount = c(a = 1, b = 2, c = 3))

a 1
b 2
c 3

Such positions can now be passed as arguments current and target into function rebalance.

We want to rebalance, and we have a vector of current weights and a vector of target weights.

The following rules apply:

• new titles (i.e. with current weight zero) are bought and get their target weight
• titles that have a zero target weight are removed from the portfolio
• after the first two rules have been applied, there will probably remain a non-zero cash position. We try to reduce it to zero with the least number of trades: when buying, we start with the assets with the lowest current weights, and vice versa

The function rebalance1 takes as input the current portfolio (\texttt{current}), the target portfolio (\texttt{target}) and min max weights (\texttt{wmin} und \texttt{wmax}). The target portfolio \texttt{target} must conform with the weight limits, i.e. all weights must be between \texttt{wmin} und \texttt{wmax}.

rebalance1 <- function(current, target, wmin = 0.025, wmax = 0.075) {
    stopifnot(length(current) == length(target))
    stopifnot(wmax >= wmin)
    zero <- 1e-10

    ans <- numeric(length(current))

    ## new assets
    i <- current < zero & target > zero
    ans[i] <- target[i]

    # old and new assets
    i <- current > zero & target > zero
    ans[i] <- pmin(pmax(current[i], wmin), wmax)

    cash <- 1 - sum(ans)

    pos <- target > zero
    while (cash > 0) {
        room <- wmax - ans[pos]
        i <- which.max(room)[1]
        eps <- min(room[i], cash)
        ans[pos][i] <- ans[pos][i] + eps
        cash <- cash - eps
    }
    while (cash < 0) {
        room <- ans[pos] - wmin
        i <- which.max(room)[1]
        eps <- min(room[i], -cash)
        ans[pos][i] <- ans[pos][i] - eps
        cash <- cash + eps
    }
    ans
}

A test: random portfolios random_p

random_p <- function(n, wmin = 0.01, wmax = 0.09) {
    .min <- 0
    .max <- 2
    while(.min < wmin || .max > wmax) {
        k <- sample(18:25,1)
        ans <- numeric(n)
        ans[sample(n,k)] <- runif(k)
        ans <- ans/sum(ans)
        .min <- min(ans[ans > 0])
        .max <- max(ans)
    }
    ans
}

current <- random_p(30)
target <- random_p(30)
new <- rebalance1(current, target)

data.frame(
    current = current,
    target = target,
    new = new,
    weights_differ = current != target,
    do_trade = current != new)

Whenever you need to round positions, you may prefer to do an actual optimisation. The ideal place for this optimisation is the original objective function, not in rebalance. And the differences, if there are any at all, are typically small. But here is an example.

n <- 10
target <- runif(n)
target <- target/sum(target)
price <- sample(10:200, n, replace = TRUE)
s <- sample(c(1,5,10,100), n, replace = TRUE,
            prob = c(0.4,0.4,0.1,0.1))
data.frame(price = price, lot.size = s)

   price lot.size
1    178        5
2     37        5
3     62        5
4     93        1
5     81        5
6    111        5
7    146        5
8    154        5
9    187        1
10   138        1

Now suppose we have only a limited budget available.

budget <- 10000
x <- rebalance(0, target, notional = budget,
               price = price, match.names = FALSE)
x

   price current value   %       target value    %       order
1    178       0     0 0.0            4   712  7.1           4
2     37       0     0 0.0           40  1480 14.8          40
3     62       0     0 0.0           20  1240 12.4          20
4     93       0     0 0.0           16  1488 14.9          16
5     81       0     0 0.0           13  1053 10.5          13
6    111       0     0 0.0            6   666  6.7           6
7    146       0     0 0.0            4   584  5.8           4
8    154       0     0 0.0            6   924  9.2           6
9    187       0     0 0.0            5   935  9.3           5
10   138       0     0 0.0            7   966  9.7           7

Notional: 10000.  Amount invested: 10048.  Total (2-way) turnover: 10048.

Now we use TAopt, from the NMOF package, to find the optimal integer portfolio.

require("NMOF")
ediff <- function(x) {
    tmp <- x*price/budget - target
    sum(tmp*tmp)
}

neighbour <- function(x) {
    i <- sample.int(length(x), size = 1L)
    x[i] <- x[i] + if (runif(1) > 0.5) - s[i] else s[i]
    x
}

sol <- TAopt(ediff,
             algo = list(x0 = numeric(length(price)),
                         neighbour = neighbour,
                         q = 0.1,
                         nS = 1000,
                         printBar = FALSE))

Threshold Accepting.

Computing thresholds ... OK.
Estimated remaining running time: 0.23 secs.

Running Threshold Accepting...
Initial solution:  0.109341
Finished.
Best solution overall: 0.001108741

df <- data.frame(TA = sol$xbest, rounded = s*round(x$target/s))
df[apply(df, 1, function(i) any(i != 0)), ]

   TA rounded
1   5       5
2  40      40
3  20      20
4  16      16
5  15      15
6   5       5
7   5       5
8   5       5
9   5       5
10  7       7

The difference.

ediff(sol$xbest) - ediff(s*round(x$target/s))

[1] 0

If you run tests with baskets of instruments or whole strategies, you often need to substitute the components of the basket for overall basket. PMwR provides a function replace_weight that helps with this task. (It is also helpful if you have hierarchies of benchmarks or want to do a `lookthrough' through a subportfolio within your portfolio.)

Suppose we have this weight vector:

w <- c(basket_1 = 0.3,
       basket_2 = 0.5,
       basket_3 = 0.2)

We also know what the first two baskets represent.

b1 <- c(a = 0.5, b = 0.2, c = 0.3)
b2 <- c(d = 0.1, e = 0.2, a = 0.7)

Now we can call replace_weight.

replace_weight(w,
               basket_1 = b1,
               basket_2 = b2)

basket_1::a basket_1::b basket_1::c
       0.15        0.06        0.09

basket_2::d basket_2::e basket_2::a
       0.05        0.10        0.35

   basket_3
       0.20

If the names of the baskets or of the things in the baskets have spaces or other characters that cause trouble, quote them.

replace_weight(c("basket 1" = 0.3,
                 "basket 2" = 0.7),
               "basket 1" = b1,
               "basket 2" = b2)

PMwR comes with a dataset called DAX, which stands for Deutscher Aktienindex (German Equity Index). The dataset is a data-frame of one column that contains the price for the day, with the timestamps stored as rownames in format YYYY-MM-DD.

str(DAX)

'data.frame': 505 obs. of  1 variable:
 $ DAX: num  9400 9435 9428 9506 9498 ...

head(DAX)

               DAX
2014-01-02 9400.04
2014-01-03 9435.15
2014-01-06 9428.00
2014-01-07 9506.20
2014-01-08 9497.84
2014-01-09 9421.61

We first transform the dataframe into an NAVseries by calling the function of the same name.

dax <- NAVseries(DAX[[1]],
                 timestamp = as.Date(row.names(DAX)),
                 title = "DAX")
dax

A concise summary is printed.

DAX
02 Jan 2014 ==> 30 Dec 2015   (505 data points, 0 NAs)
    9400.04           10743

There is also a generic function as.NAVseries, which can be used, for instance, to coerce zoo series to NAVseries. Or, after having run a backtest (see Chapter Backtesting), saying

as.NAVseries(btest(....))

extracts the NAV series from the backtest data.

Calling summary on an NAV series should provide useful statistics of the series. You may notice that the summary does not provide too many statistics. It has been said that the purpose of statistics is to reduce many numbers to only few – not the other way around. Thus, summary.NAVseries will confine itself to few statistics that can be computed in a reasonably robust way. By that I mean that a statistic is useful for many types of NAV series.¹¹ Numbers that are provided should be formatted and presented in a way that is inline with (good) industry practice and makes sense statistically. Returns, for instance, will only be annualised when the NAV series spans more than one calendar year; volatility will always be annualised and be computed from monthly data (if possible); all numbers are rounded to a meaningful precision \citep{Ehrenberg1981}.

summary(dax)

---------------------------------------------------------
DAX
02 Jan 2014 ==> 30 Dec 2015  (505 data points, 0 NAs)
    9400.04           10743
---------------------------------------------------------
High               12374.73  (10 Apr 2015)
Low                 8571.95  (15 Oct 2014)
---------------------------------------------------------
Return (%)              6.9  (annualised)
---------------------------------------------------------
Max. drawdown (%)      23.8
_ peak             12374.73  (10 Apr 2015)
_ trough            9427.64  (24 Sep 2015)
_ recovery                   (NA)
_ underwater now (%)   13.2
---------------------------------------------------------
Volatility (%)         18.0  (annualised)
_ upside               14.4
_ downside             10.4
---------------------------------------------------------

Monthly returns  ▁▆▇▃▆▃▁

      Jan Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct Nov  Dec YTD
2014 -1.0 4.1 -1.4  0.5  3.5 -1.1 -4.3  0.7  0.0 -1.6 7.0 -1.8 4.3
2015  9.1 6.6  5.0 -4.3 -0.4 -4.1  3.3 -9.3 -5.8 12.3 4.9 -5.6 9.6

For summaries of NAV series, a method for toLatex can be used to fill \LaTeX-templates. The package comes with a vignette that provides examples.

PMwR offers two functions that may provide insights into NAV series: drawdowns and streaks. A vignette shows examples.

vignette("Drawdowns_streaks", package = "PMwR")

To explain what the function does, we use two very short time-series: the values of the DAX, the German stock-market index, and the REXP, a German government-bond index, from 2 January and 8 January 2014 (just 5 trading days). We also combine them into a matrix drax.

dax  <-  DAX[1:5, ]
rexp <- REXP[1:5, ]
drax <- cbind(dax, rexp)

Calling scale1 on dax is equivalent to dividing the whole series by its first element.

scale1(dax) ==  dax/dax[1]

[1] TRUE TRUE TRUE TRUE TRUE

Lest you skip the rest of the chapter: scale1 comes with several additional features.

It is common, too, to scale to a level of 100. We either multiply the whole series by 100, or use the level argument.

scale1(dax, level = 100)

[1] 100.00 100.37 100.30 101.13 101.04
attr(,"scale1_origin")
[1] 1

(The scale1_origin attribute will be explained shortly.)

If we give a matrix to scale1, the function scales each column separately.

scale1(drax, level = 100)

        dax   rexp
[1,] 100.00 100.00
[2,] 100.37 100.06
[3,] 100.30 100.23
[4,] 101.13 100.29
[5,] 101.04 100.28
attr(,"scale1_origin")
[1] 1

scale1 is a generic function; it works, for instance, with zoo objects.

library("zoo")
drax.zoo <- zoo(drax, as.Date(row.names(DAX)[1:5]))
scale1(drax.zoo, level = 100)

              dax   rexp
2014-01-02 100.00 100.00
2014-01-03 100.37 100.06
2014-01-06 100.30 100.23
2014-01-07 101.13 100.29
2014-01-08 101.04 100.28
attr(,"scale1_origin")
[1] 2014-01-02

plot(scale1(drax.zoo, level = 100),
     plot.type = "single",
     xlab = "",
     ylab = "",
     col = c("darkblue", "darkgreen"))

The argument when defines the origin.

scale1(drax, when = 3, level = 100)

         dax    rexp
[1,]  99.703  99.769
[2,] 100.076  99.830
[3,] 100.000 100.000
[4,] 100.829 100.062
[5,] 100.741 100.049
attr(,"scale1_origin")
[1] 3

This origin is attached to the scaled series as an attribute scale1_origin. This is useful if you want mark the start of the scaled series; for instance, in a plot with abline.

With a zoo object, when should be compatible with the class of the object's index.

scale1(drax.zoo, when = as.Date("2014-01-07"), level = 100)

               dax    rexp
2014-01-02  98.883  99.707
2014-01-03  99.253  99.768
2014-01-06  99.177  99.938
2014-01-07 100.000 100.000
2014-01-08  99.912  99.987
attr(,"scale1_origin")
[1] 2014-01-07

when also understands the keyword first.complete, which is actually the default, and the keywords first and last. first.complete is useful when some series have leading missing values.

drax[1:2, 1] <- NA
drax

        dax   rexp
[1,]     NA 440.53
[2,]     NA 440.79
[3,] 9428.0 441.55
[4,] 9506.2 441.82
[5,] 9497.8 441.76

scale1(drax, level = 100)   ## 'first.complete' is the default

        dax    rexp
[1,]     NA  99.769
[2,]     NA  99.830
[3,] 100.00 100.000
[4,] 100.83 100.062
[5,] 100.74 100.049
attr(,"scale1_origin")
[1] 3

When the argument centre is TRUE, the mean return is subtracted from the returns.

scale1(drax.zoo, centre = TRUE)

               dax    rexp
2014-01-02 1.00000 1.00000
2014-01-03 1.00114 0.99991
2014-01-06 0.99779 1.00091
2014-01-07 1.00348 1.00083
2014-01-08 1.00000 1.00000
attr(,"scale1_origin")
[1] 2014-01-02

The default is to subtract the geometric mean: the series will have a growth rate of zero; it will end where it started.

The argument scale takes a standard deviation and scales the returns to that standard deviation.

apply(returns(scale1(drax.zoo, scale = 0.02)), 2, sd)

 dax rexp
0.02 0.02

This may create fairer comparisons, for instance, between fund prices that exhibit very different volatilities. It can also help to visualise correlation.

It should be stressed that centre and scale change returns, but scale1 expects and evaluates to levels (not returns).

The zoo method has a further argument that affects returns: inflate. To illustrate its use, let us create a constant series.

z <- zoo(100,
         seq(from = as.Date("2015-1-1"),
             to   = as.Date("2016-1-1"),
             by   = "1 day"))
head(z)
tail(z)

2015-01-01 2015-01-02 2015-01-03 2015-01-04 2015-01-05 2015-01-06
       100        100        100        100        100        100

2015-12-27 2015-12-28 2015-12-29 2015-12-30 2015-12-31 2016-01-01
       100        100        100        100        100        100

inflate should be a numeric value: the annual growth rate that is added to the time-series's return (or that is subtracted from it, if negative).

head(scale1(z, inflate = 0.02))
tail(scale1(z, inflate = 0.02))

2015-01-01 2015-01-02 2015-01-03 2015-01-04 2015-01-05 2015-01-06
    1.0000     1.0001     1.0001     1.0002     1.0002     1.0003

2015-12-27 2015-12-28 2015-12-29 2015-12-30 2015-12-31 2016-01-01
    1.0197     1.0198     1.0198     1.0199     1.0199     1.0200

The previous section provided examples of scaling series. In this section, we are going to see how scale1 does its computations.

First, a series \(P\) passed to scale1 is transformed into returns, \(R\). The scale argument allows you to set a desired volatility for the series's returns, defined as their standard deviation. The computation uses the fact that multiplying a random variable by a number \(b\) changes its variance to \(b^2\) times its original variance. Hence, scale1 divides the returns by the actual standard deviation and then multiplies them by the desired one (i.e. the value passed via the scale argument).

Changing total return (or, equivalently, average return) is slightly more complicated. Suppose we want to scale the total return of the series \(P\) such that it equals some target return \(r_*\). Start with writing the total return as the product of single-period returns.

There clearly is an infinity of possible adjustments that would do the trick. We might, for instance, change only \(P_0\) or \(P_T\) so that the desired return is achieved.

But that is probably not what we want. A reasonable requirement is that the scaling touches as few other statistical properties as possible. Adding a constant \(z\) to the return in every period does that: it does not change the volatility of the returns; neither does it affect linear or rank correlation of the returns with some other series. To compute \(z\), we need to solve the following equation:

Alternatively, we may use logs.

This is an application for root-finding (see chapter 11 of \citealp{Gilli2019}), for which scale1 uses uniroot.

We have the following trades and times.

amount <- c(1, 3, -3, 1, -3,  1)
time   <- c(0, 1,  3, 4,  7, 12)

The holding period (duration) of these trades can be computed so:

data.frame(position = cumsum(amount)[-length(amount)],
           from = time[-length(time)],
           to   = time[-1L],
           duration = diff(time))

  position from to duration
1        1    0  1        1
2        4    1  3        2
3        1    3  4        1
4        2    4  7        3
5       -1    7 12        5

We can plot the exposure.

plot(c(time[1], time), cumsum(c(0, amount)),
     type = "s", xlab = "time", ylab = "position")

Thus, we have had a position from time 0 to time 12 (hours into the trading day, say), but its size varied. The function tw_exposure (time-weighted exposure) computes the average absolute exposure.

tw_exposure(amount, time)

1.75

To give a simpler example: suppose we bought at the open of a trading day and sold at noon. The average exposure for the day is thus half a contract.

amount <- c(1,  -1, 0)
time   <- c(0, 0.5, 1)
tw_exposure(amount, time)

0.5

If we bought at the open, went short at noon, and closed the position at the end of the day, the average exposure would be one contract, since absolute position size is relevant.

amount <- c(1, -2 , 1)
time   <- c(0,0.5,1)
tw_exposure(amount, time)

1

Whether absolute exposure is used is controlled by an argument abs.value. Setting it to FALSE can be useful to detect long or short biases.

tw_exposure(amount, time, abs.value = FALSE)

0

We have the following trades.

timestamp <- 1:3
amount <- c(-1, 2, -1)
price <- c(100, 99, 101)

Calling split_trades will return a list of two single trades. Each single trade, in turn, is a list with components amount, price and timestamp.

split_trades(amount = amount,
             price = price,
             timestamp = timestamp,
             aggregate = FALSE)

[[1]]
[[1]]$amount
[1] -1  1

[[1]]$price
[1] 100  99

[[1]]$timestamp
[1] 1 2


[[2]]
[[2]]$amount
[1]  1 -1

[[2]]$price
[1]  99 101

[[2]]$timestamp
[1] 2 3

Note that the second transaction (buy 2 @ 99) has been split up: buying one contract closes the first trade; the other contract opens the second trade. This splitting is useful in its own right: there are accounting systems around that cannot handle a trade that switches a position directly from long to short, or vice versa. Instead, the trade needs first be closed (i.e. the net position becomes zero).

With argument aggregate set to TRUE, the function reconstructs the total series, but with those trades splitted that change the position's sign.

split_trades(amount, price, timestamp, aggregate = TRUE)

$amount
[1] -1  1  1 -1

$price
[1] 100  99  99 101

$timestamp
[1] 1 2 2 3

Another example. We have the following trades and impose a limit that the maximum absolute exposure for the trader should only be 2.

timestamp <- 1:6
amount <- c(-1,-1,-1,1,1,1)
price <- c(100,99,98,98,99,100)
limit(amount, price, timestamp, lim = 2)

$amount
[1] -1 -1  1  1

$price
[1] 100  99  99 100

$timestamp
[1] 1 2 5 6

Scaling the trades.

scale_to_unity(amount)

[1] -0.333 -0.333 -0.333  0.333  0.333  0.333

Closing the trade at once.

close_on_first(amount)

[1] -1 -1 -1  3  0  0

We have the following sample of prices of the Bund future contract, traded at the Eurex in Germany.

(You'll find the code to generate those times in the code file for this chapter.)

Note that I have left the time zone to the operating system. Since my computer is typically located in the time zone that the tz database (http://www.iana.org/time-zones)

calls 'Europe/Zurich', the first time should be 2012-10-18 20:00:09. If, for instance, your computer is in 'America/Chicago' instead and you run the above code, the first time would be 2012-10-18 13:00:09. Which is right: it is the correct time, only translated into Chicago local time.

A plot of price against time looks like this.

plot(times, prices, type = "s")

Such a plot is fine for many purposes. But the contract for which we have prices is only traded from Monday to Friday, not on weekends, and it is traded only from 08:00 to 22:00 Europe/Berlin time. So the plot should omit those times at which no trading takes place. This is what the function plot_trading_hours does.

tmp <- plot_trading_hours(x = prices, t = times,
                          interval = "1 sec",
                          labels = "day",
                          fromHHMMSS = "080000",
                          toHHMMSS = "220000",
                          type = "s")

What we need for such a plot is a function that maps actual time to a point on the x-scale, while the y-scale stays unchanged. If we were talking only about days, not times, we needed something like this:

This mapping is what plot_trading_hours provides. And not much more: the design goal of the function is to make it as much as possible an ordinary plot; or more specifically, to make it as similar as possible to the plot function. Indeed, plot_trading_hours calls plot with a small number of default settings:

list(type = "l", xaxt = "n", xlab = "", ylab = "")

These settings can all be overridden through the ... argument, which is passed to plot. Note that we already set s as the plot's type in the last code chunk. The only required setting is suppressing the x-axis with setting xaxt to 'n', because plot_trading_hours will create its own x-axis via a call to axis(1, ...). In case you wish to use your own axis specification, either set do.plotAxis to FALSE or pass settings to axis through the list axis1.par.

plot_trading_hours also handles data on a daily frequency. The function will assume such a case if the timestamp is of class Date; it will then choose sensible defaults for the time-axis. In effect, the function will remove weekends and, if specified, holidays.

As an example, consider the first 10 observations of the DAX dataset.

x <- DAX[1:10, ]
t <- as.Date(row.names(DAX)[1:10])
data.frame(t, x, weekday = weekdays(t))

If we plot these data, there will be a value plotted for 2014-01-04, as marked by the vertical line.

plot(t, x, type = "l")
abline(v = as.Date("2014-01-04"))

That is despite the fact that this is a Saturday.

format(as.Date("2014-01-04"), "%A")

Saturday

This is not a bug: it is the default behaviour of plot. Saturday's value results from an interpolation. plot_trading_hours instead ignores such nonexisting dates.

plot_trading_hours(x, t)

Valuing an instrument can mean using either a market price or a theoretical price. In the discussion that follows, I will assume that we already have prices (or net-present values).

The function div_adjust corrects price series for dividends. It is meant as a low-level function and is implemented to work on numeric vectors. Consider a hypothetical price series x, which goes ex-dividend at time 3.

x <- c(9.777, 10.04, 9.207, 9.406)
div <- 0.7
t <- 3

The default for div_adjust is to match the final price.

div_adjust(x, t, div)

[1] 9.086185 9.330603 9.207000 9.406000

If you prefer a correction that matches the first price, set argument backward to FALSE.

div_adjust(x, t, div, backward = FALSE)

[1]  9.77700 10.04000  9.90700 10.12113

The function split_adjust handles stock splits. It is implemented to work on numeric vectors.

US treasury bonds are often quoted in 1/32nds of points. For instance, the price 110'030 would mean 110+3/32.

The function quote32 lets you parse and `pretty-print' such prices.

quote32(c("110-235", "110-237"))

[1] 110-23+ 110-23¾

Internally, quote32 will store the prices as numeric values: the fractions are only used for printing.

as.numeric(quote32(c("110-235", "110-237")))

[1] 110.73 110.74

dput(quote32(c("110-235", "110-237")))

structure(c(110.734375, 110.7421875),
          handle = c(110, 110),
          ticks = c(23, 23),
          fraction = c(2, 3),
          class = "quote32")

An ISIN, which stands for International Securities Identification Number, uniquely¹³ identifies a security.

is_valid_ISIN(c("DE0007236101",     ## Siemens
                "DE0007236102"))    ## last digit changed

[1]  TRUE FALSE

There is a function is\_valid\_SEDOL too.

A pricetable is a matrix of prices, with some added functionality for subsetting.

To normal people, a tree consists of a trunc, branches and leaves. To people who do graph theory, a tree is a connected graph with only one path between any two nodes.

Trees are useful to represent hierachies – just think of a file tree. For portfolios, a tree can be used to indicate groupings, such as countries or industries. (See argument account in function position.)

Methods are responsible for 'stripping' the input down do x and t, calling 'returns.default' or some other method, and then to re-assemble the original class's structure. When period is not specified, methods should keep timestamp information for themselves and not pass it on. That is, returns.default should only ever receive a timestamp when period is specified.