Trying to do a simulation in R

Question

I'm pretty new to R, so I hope you can help me! I'm trying to do a simulation for my Bachelor's thesis, where I want to simulate how a stock evolves. I've done the simulation in Excel, but the problem is that I can't make that large of a simulation, as the program crashes! Therefore I'm trying in R.

The stock evolves as follows (everything except $\epsilon$ consists of constants which are known): $$W_{t+\Delta t} = W_t exp^{r \Delta t}(1+\pi(exp((\sigma \lambda -0.5\sigma^2) \Delta t+\sigma \epsilon_{t+\Delta t} \sqrt{\Delta t}-1))$$

The only thing here which is stochastic is $\epsilon$, which is represented by a Brownian motion with N(0,1).

What I've done in Excel:

• Made 100 samples with a size of 40. All these samples are standard normal distributed: N(0,1).
• Then these outcomes are used to calculate how the stock is affected from these (the normal distribution represent the shocks from the economy).

My problem in R:

I've used the sample function: x <- sample(norm(0,1), 1000, T) So I have 1000 samples, which are normally distributed. Now I don't know how to put these results into the formula I have for the evolution of my stock. Can anyone help?

Answer 1

Using R for (discrete) simulation

There are two aspects to your question: conceptual and coding.

Let's deal with the conceptual first, starting with the meaning of your equation:

1. Conceptual issues

The first thing to note is that your evolution equation is continuous in time, so running your simulation as described above means accepting a discretisation of the problem. Whether or not that is appropriate depends on your model and how you have obtained the evolution equation.

If you do run a discrete simulation, then the key decision you have to make is what stepsize $\Delta t$ you will use. You can explore different step-sizes to observe the effect of step-size, or you can proceed analytically and attempt to derive an appropriate step-size.

Once you have your step-size, your simulation consists of pulling new shocks (samples of your standard normal distribution), and evolving the equation iteratively until the desired time has elapsed. The final state $W_t$ is then available for you to analyse however you wish. (If you retain all of the $W_t$, you have a distribution of the trajectory of the system as well, which you can analyse.)

So:

• your $x$ are a sampled distribution of your shocks, i.e. they are $\epsilon_t=0$.
• To simulate the evolution of the $W_t$, you will need some initial condition $W_0$. What this is depends on what you're modelling. If you're modelling the likely values of a single stock starting at an initial price $W_0$, then your initial state is a 1000 element vector with constant value.
• Now evaluate your equation, plugging in all your constants, $W_0$, and your initial shocks $\epsilon_0 = x$ to get the distribution of prices $W_1$.
• Repeat: sample $x$ again -- this is now $\epsilon_1$. Plugging this in, gives you $W_2$ etc.

2. Coding the simulation (simple example)

One of the useful features of R is that most operators work element-wise over vectors.

So you can pretty much type in your equation more or less as it is.

I've made a few assumptions about the parameters in your equation, and I've ignored the $\pi$ function -- you can add that in later.

So you end up with code that looks something like this:

 dt <- 0.5   # step-size
 r <- 1      # parameters
 lambda <- 1
 sigma <- 1    # std deviation
 w0 <- rep(1,1000)   # presumed initial condition -- prices start at 1

 # Show an example iteration -- incorporate into one line for production code...
 x <- rnorm(1000,mean=0,sd=1)  # random shock
 w1 <- w0*exp(r*dt)*(1+exp((sigma*lambda-0.5*sigma^2)*dt +
           sigma*x*sqrt(dt) -1))  # evolution

When you're ready to let the simulation run, then merge the last two lines, i.e. include the sampling statement in the evolution statement. You then get one line of code which you can run manually or embed into a loop, along with any other analysis you want to run.

 # General simulation step
 w <- w*exp(r*dt)*(1+exp((sigma*lambda-0.5*sigma^2)*dt +  
           sigma*rnorm(1000,mean=0,sd=1)*sqrt(dt) -1))

You can also easily visualise the changes and obtain summary statistics (5-number summary):

 hist(w)
 summary(w)

Of course, you'll still need to work through the details of what you actually want to model and how you want to go about analysing it --- and you've got the $\pi$ function to deal with --- but this should get you started toward using R for discrete simulation.