#install.packages("pbapply")
library(quantmod)
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: TTR
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
#install.packages("lubridate")
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(xts)
library(zoo)
library(pbapply)

# Define the ticker symbol and download data
getSymbols("AMZN", from = "2000-01-01", to = Sys.Date())
## [1] "AMZN"
# Use Ad() to get Adjusted closing prices
tmp <- Ad(AMZN)

# Print the adjusted prices
print(tmp)
##            AMZN.Adjusted
## 2000-01-03      4.468750
## 2000-01-04      4.096875
## 2000-01-05      3.487500
## 2000-01-06      3.278125
## 2000-01-07      3.478125
## 2000-01-10      3.459375
## 2000-01-11      3.337500
## 2000-01-12      3.178125
## 2000-01-13      3.296875
## 2000-01-14      3.212500
##        ...              
## 2024-02-20    167.080002
## 2024-02-21    168.589996
## 2024-02-22    174.580002
## 2024-02-23    174.990005
## 2024-02-26    174.729996
## 2024-02-27    173.539993
## 2024-02-28    173.160004
## 2024-02-29    176.759995
## 2024-03-01    178.220001
## 2024-03-04    177.580002
# Create a variable to store the returns
rets <- ROC(tmp, type = "discrete")
rets[is.na(rets)] <-0
mean(rets)
## [1] 0.001092086
sd(rets)
## [1] 0.03145133
# Create a formula to return the predicted price
stk_ret = function(STK_PRC, N, MEAN, STDEV)
{
  delta_t = 1/N # for 1 period
  for (i in seq(N))
  {
    epsilon <- runif(n=1, min=0, max=1) # random probabilities
    STK_PRC <- STK_PRC * (1+ qnorm(epsilon, MEAN*delta_t, STDEV*sqrt(delta_t)))
  }
  STK_PRC
}

#latest quote
last(tmp)
##            AMZN.Adjusted
## 2024-03-04        177.58
library(lubridate)

#run 1000 simulations
simulations <- 1000

# number of trading days
N = 20 #look back 20 days to predict the stock price, this will be our window

STK_PRC <- as.numeric(coredata(tmp[Sys.Date() - days(20)]))

# Get the mean and the standard deviation of the returns
MEAN = mean(rets)
STDEV = sd(rets)

#Print the results
MEAN
## [1] 0.001092086
STDEV
## [1] 0.03145133
# Create a vector to store stock prices
stock_prices <- c()
for (i in seq(simulations))
{
  stock_prices <- c(stock_prices, stk_ret(STK_PRC = STK_PRC, N=N, MEAN=MEAN, STDEV = STDEV))
}

# Show the predictions
stock_prices
##    [1] 179.8883 173.3379 167.5631 178.1745 172.1231 169.2847 168.4716 171.7750
##    [9] 172.6127 167.8123 168.7825 166.4053 166.3646 187.0783 169.9612 174.8407
##   [17] 166.1696 182.4276 173.5382 167.5846 176.1459 169.7327 169.0370 174.1323
##   [25] 173.5868 172.8030 168.6842 165.7751 165.5367 169.1753 168.3848 167.0489
##   [33] 173.4539 170.5004 167.1037 169.4842 168.4603 175.5009 172.2202 174.7537
##   [41] 174.0878 174.3028 166.0618 170.7904 180.0490 170.2998 175.7720 170.4757
##   [49] 173.0875 173.0097 174.5601 183.9669 170.4239 169.6498 173.0566 166.3239
##   [57] 166.1644 164.3771 168.4704 173.7635 161.1323 175.3631 165.6393 171.4321
##   [65] 168.9418 178.4871 163.7741 171.5211 165.4294 169.8396 173.8939 178.4756
##   [73] 176.9219 173.1431 165.5694 168.4637 178.2430 169.3285 179.1799 161.2057
##   [81] 171.3836 169.1949 175.3589 164.5995 154.5033 177.6439 164.0728 174.7456
##   [89] 176.8971 164.3284 175.3681 173.7958 166.6729 175.4211 181.2745 168.0598
##   [97] 171.0741 174.7974 157.4639 181.1324 181.8000 181.3308 171.4393 171.7618
##  [105] 173.8863 169.0829 180.1170 172.1168 175.2433 167.0541 163.0638 167.4174
##  [113] 170.6878 177.4565 168.6547 169.9070 174.6807 171.3280 170.2063 174.2719
##  [121] 179.5230 169.8973 166.7248 167.1808 170.9291 169.3130 175.3940 172.1737
##  [129] 155.0036 167.6781 160.8006 169.9078 162.0489 166.0214 183.4650 175.3754
##  [137] 174.9630 176.6373 177.7469 170.8460 165.1088 174.7537 178.0822 170.2971
##  [145] 172.8976 169.8177 179.2994 170.7951 178.0476 167.7336 174.5495 163.9813
##  [153] 174.3910 163.9054 165.3527 163.3891 160.0488 164.3234 174.9069 173.3275
##  [161] 166.5272 174.7909 175.0597 167.2875 164.6784 172.4353 175.3906 175.3698
##  [169] 168.0654 166.5359 176.8767 173.1856 178.0717 165.0677 174.6684 170.2893
##  [177] 174.0220 174.4307 167.1545 172.1824 174.2902 174.2278 173.5633 167.1189
##  [185] 162.8769 179.1087 166.2088 178.7012 161.9299 165.6654 181.2937 168.8133
##  [193] 161.0979 172.3312 174.0108 170.3525 160.9274 179.5956 163.5928 175.1496
##  [201] 166.3548 170.5923 178.2079 165.8085 168.6613 173.9860 170.6243 166.4569
##  [209] 166.4505 177.2566 173.3842 176.1825 170.1342 168.7947 179.6524 175.1227
##  [217] 172.2294 173.8626 175.0727 171.6023 168.6695 177.2964 164.3908 167.1910
##  [225] 180.2779 162.0318 167.8611 164.6187 167.5849 166.2809 161.5713 169.2913
##  [233] 173.5382 173.3303 175.5111 178.9047 176.1741 165.0819 171.7976 171.7231
##  [241] 180.6626 170.2193 171.1419 166.8245 170.4848 168.1022 174.4102 172.6890
##  [249] 166.7978 168.6434 166.8854 164.8882 164.1926 172.9631 171.4603 174.1183
##  [257] 168.0783 171.7180 167.6075 171.5463 178.1676 167.5134 175.5863 174.1372
##  [265] 174.2056 170.6442 162.8383 169.4766 160.9505 181.4141 163.6525 177.1607
##  [273] 172.8172 171.7948 165.9260 171.1283 171.3052 179.8329 169.9239 162.8864
##  [281] 175.0224 175.1508 169.9282 184.6126 164.2549 171.8451 170.4118 177.8098
##  [289] 171.7822 165.7574 184.4672 176.3061 174.7943 157.0441 172.0286 180.3842
##  [297] 174.9488 175.1283 179.2561 173.9041 173.1053 179.3675 178.3865 162.9654
##  [305] 166.4697 170.3541 166.8783 170.8628 174.6360 164.5319 169.6786 180.4904
##  [313] 175.3941 170.7595 172.4605 166.9031 163.5994 175.0448 158.9443 178.2024
##  [321] 173.4686 162.9623 176.3859 170.6236 175.4447 176.2052 171.1974 168.0246
##  [329] 179.0568 186.7956 172.7640 175.6441 170.1568 169.1344 166.7790 174.3026
##  [337] 175.3603 164.1287 171.2801 167.1170 174.6383 169.8601 175.2558 164.0798
##  [345] 168.2031 164.9812 175.8983 173.6563 164.1802 171.9243 171.6233 168.2180
##  [353] 175.1815 166.7971 167.8582 167.3244 165.7315 185.4097 168.9546 172.1459
##  [361] 175.6734 167.0686 162.3954 159.9886 167.1188 174.4437 187.4702 167.2950
##  [369] 166.9357 164.7651 171.6023 180.9847 174.9420 165.6122 175.2689 166.6059
##  [377] 174.0676 169.5048 170.3016 180.0108 167.8215 178.0424 176.5589 161.8344
##  [385] 174.2813 172.3375 175.0179 174.0617 177.6385 173.9916 178.6602 169.6451
##  [393] 177.0683 170.0849 178.5939 164.4373 173.0822 166.6530 165.9618 170.4366
##  [401] 174.9498 168.7132 163.8575 182.2717 168.2649 170.6847 170.1502 168.0557
##  [409] 162.4432 174.8292 168.3537 165.6196 175.8656 172.7026 176.7250 170.0864
##  [417] 170.6313 171.4338 172.9179 170.6131 173.5479 178.7912 170.4670 165.3618
##  [425] 178.6768 174.2956 163.8295 172.3219 182.2328 170.1424 163.0891 156.8658
##  [433] 168.5879 167.0860 175.3183 165.6323 173.7643 173.0985 168.9281 166.0713
##  [441] 165.6342 185.5429 177.2794 173.3030 169.4267 174.4487 170.6779 178.6884
##  [449] 168.3155 173.4165 168.7744 172.2947 173.3295 161.6078 173.1651 163.8758
##  [457] 171.6968 177.5830 179.3342 165.2755 164.2780 171.4702 169.8845 168.8325
##  [465] 177.2288 167.6772 167.0125 172.5998 173.5111 167.8089 173.2020 168.8035
##  [473] 171.5827 165.1086 174.4130 171.8647 170.6986 175.6320 176.8819 168.5918
##  [481] 172.0383 176.7991 177.5071 178.7243 161.9400 171.1487 176.9460 164.4290
##  [489] 171.0103 175.7048 178.7481 178.1342 172.1676 166.7504 173.9783 169.0999
##  [497] 172.8192 172.8870 178.0868 178.1846 172.6552 179.8461 181.8988 180.9349
##  [505] 163.6540 174.8240 161.5459 164.5370 168.5552 163.4750 159.8253 168.6006
##  [513] 169.1458 169.2838 167.9351 187.3650 166.9259 178.1290 162.6400 171.6135
##  [521] 164.7202 167.0336 169.7884 171.1588 170.3439 174.1473 175.1174 177.6526
##  [529] 170.0033 169.1282 180.0809 177.0640 162.2770 163.2760 165.3118 172.3621
##  [537] 168.2384 172.3912 175.7491 174.8843 167.3223 162.3369 170.5510 168.9989
##  [545] 169.9185 168.3935 169.7825 172.4269 164.2683 175.9633 167.5285 164.4618
##  [553] 166.8056 179.1277 168.7767 168.7322 182.7701 171.0126 165.9911 170.3856
##  [561] 165.2762 170.5153 172.9208 174.7160 176.4739 175.4325 162.8680 170.3109
##  [569] 176.6642 172.8222 168.3471 172.8632 172.2634 173.2023 175.0967 176.2518
##  [577] 167.3326 164.2942 169.2280 167.6074 179.3212 167.2639 173.6614 174.3740
##  [585] 167.3783 169.4811 166.8409 164.6919 177.4955 167.0919 172.7442 170.0497
##  [593] 169.8879 166.6744 177.7046 168.4651 172.8293 167.2691 172.9239 180.2325
##  [601] 173.7441 173.6037 167.6144 176.4538 188.8284 169.6253 168.7328 173.8770
##  [609] 174.9406 173.7216 168.6297 175.8220 170.1954 166.6901 164.4381 164.0921
##  [617] 183.5113 178.7078 173.3990 163.1957 171.1903 174.2939 168.8603 178.8268
##  [625] 170.1483 175.4025 187.0250 174.6805 168.3667 168.1307 169.3756 169.6262
##  [633] 175.7301 167.2493 174.7850 182.4914 168.9317 162.0827 180.0506 174.8932
##  [641] 176.7439 167.3722 172.5128 176.6230 172.2206 167.1223 167.2336 160.5744
##  [649] 172.3081 172.4067 170.1169 166.6099 168.5687 169.3570 171.1263 163.9855
##  [657] 169.2790 172.5752 178.1735 172.2550 172.9049 169.5660 168.1522 181.8392
##  [665] 172.9469 175.1061 168.2041 174.1143 168.3004 172.2867 163.5966 159.8122
##  [673] 173.4384 171.0049 170.2981 177.5149 168.5197 167.1469 172.9189 181.1964
##  [681] 166.7731 163.8104 175.9785 167.2032 162.5891 173.6106 169.5697 164.6631
##  [689] 175.7653 163.9547 176.1804 169.8986 165.5538 176.1963 173.6230 180.5410
##  [697] 168.2468 175.5213 165.7722 172.2819 163.0010 173.6816 164.2559 177.1815
##  [705] 167.4945 176.8633 174.5012 182.4884 168.4792 167.2566 165.1116 160.7819
##  [713] 173.5181 172.9358 167.8014 172.0022 163.1428 166.1754 172.1477 167.3218
##  [721] 165.2814 173.7031 168.9138 171.3556 177.8646 166.1836 164.8790 173.2899
##  [729] 166.2716 168.6793 170.1483 173.7501 173.4167 164.2211 160.1836 176.2337
##  [737] 174.2246 162.3530 165.4564 168.5026 168.5075 170.6661 166.9864 163.0317
##  [745] 179.5654 178.3650 166.4345 173.6608 164.8691 166.6421 169.1666 168.1108
##  [753] 179.3340 174.3627 169.8309 163.7057 167.6293 178.7316 173.3338 160.6870
##  [761] 174.7009 179.0387 174.8926 174.2645 168.3984 172.0121 166.7618 167.4204
##  [769] 180.7992 169.2100 172.9354 171.9464 172.4451 163.8318 167.8028 164.6775
##  [777] 166.0650 171.1300 167.5848 176.2850 161.7714 171.8306 181.1648 176.4904
##  [785] 169.4122 167.9886 171.0025 168.4258 176.9000 178.9038 160.3279 178.8381
##  [793] 173.4771 173.7571 164.7509 172.8280 181.1061 178.1881 175.1581 175.6450
##  [801] 175.2296 163.9146 168.1046 168.4592 170.2506 167.4793 167.2731 169.5637
##  [809] 172.5504 169.0527 169.1160 171.3389 169.3166 163.0511 170.4058 171.3447
##  [817] 163.7499 174.4127 171.0669 166.3439 170.5027 178.4025 178.3637 167.7895
##  [825] 174.6618 171.2856 175.4895 172.7927 166.1573 172.0140 167.6985 166.9676
##  [833] 178.6358 167.6119 182.2640 177.0787 181.0174 180.8176 165.8479 170.7915
##  [841] 175.7354 174.2532 169.6666 172.8622 172.9374 175.5399 173.6226 167.7594
##  [849] 176.2455 166.1885 166.8206 170.9126 162.6184 169.1544 174.3277 170.1616
##  [857] 164.2314 167.3797 164.4949 175.6495 170.9260 175.3265 168.2671 178.2562
##  [865] 175.4075 165.5792 172.3807 168.9884 168.9882 178.2116 173.2511 160.8267
##  [873] 176.4311 168.1514 164.1734 164.1759 179.0568 171.3715 180.7116 172.3077
##  [881] 162.8155 168.5361 174.0367 172.0470 165.0025 167.8742 178.8303 169.8481
##  [889] 170.7962 166.8406 175.0261 169.8190 163.9432 174.8302 169.4296 183.0787
##  [897] 173.2909 175.6221 180.8464 166.6234 170.3090 170.1396 164.2805 177.2228
##  [905] 184.2967 168.6034 169.5785 165.0700 174.3356 173.8762 166.2639 167.4374
##  [913] 167.9913 171.6196 180.6357 166.4943 165.8385 165.5860 167.4105 180.3086
##  [921] 173.0072 176.1987 169.6341 174.6863 169.2360 172.7245 174.5700 174.3754
##  [929] 176.1703 174.0345 169.8989 176.4206 175.5010 173.0189 170.2523 176.8158
##  [937] 166.4655 171.2783 181.0085 170.2588 173.6098 165.9064 170.4065 167.8267
##  [945] 172.0920 174.2802 169.8786 169.8929 174.0067 169.9116 179.5237 168.5623
##  [953] 170.6411 163.8199 166.1164 171.7354 170.1251 169.9910 173.3435 166.1552
##  [961] 173.5294 170.2384 169.1526 175.4327 171.4299 171.6255 171.7359 168.4864
##  [969] 168.9962 173.2257 171.7686 174.5261 176.0953 177.6995 169.8226 173.3206
##  [977] 173.9910 167.1849 178.3458 177.1521 168.6152 177.3918 169.9354 170.2491
##  [985] 164.1388 165.8972 175.2886 172.9241 173.0474 160.4555 171.8847 167.6967
##  [993] 174.6547 175.6217 167.2886 166.4113 166.4758 163.2330 165.9125 176.0171
# Create a summary for the predictions
quantile(stock_prices)
##       0%      25%      50%      75%     100% 
## 154.5033 167.3793 171.0115 174.8040 188.8284

#The lowest possible value $157 and highest $192 for today’s predictions

# Get the actual values
last(tmp)
##            AMZN.Adjusted
## 2024-03-04        177.58

#Because we use the mean and the standard deviation to predict the stock prices our estimates are off. These numbers do not adjust to market trends so it does not capture the volatility and stock price movements which leads to a higher spread between the predictions and the actual stock price. The standard deviation will maintain the prices within the same range.

# Now we are going to calculate the predictions for each month

# we will use options expiration dates
EXPIRY <- tmp[options.expiry(tmp)]
EXPIRY #monthly options expiration date along with the prices
##            AMZN.Adjusted
## 2000-01-21      3.103125
## 2000-02-18      3.237500
## 2000-03-17      3.240625
## 2000-05-19      2.631250
## 2000-06-16      2.300000
## 2000-07-21      2.056250
## 2000-08-18      1.950000
## 2000-09-15      2.181250
## 2000-10-20      1.540625
## 2000-11-17      1.371875
##        ...              
## 2023-05-19    116.250000
## 2023-06-16    125.489998
## 2023-07-21    130.000000
## 2023-08-18    133.220001
## 2023-09-15    140.389999
## 2023-10-20    125.169998
## 2023-11-17    145.179993
## 2023-12-15    149.970001
## 2024-01-19    155.339996
## 2024-02-16    169.509995
# Now we will try to predict the option expiration price for each date
EXPIRY <- EXPIRY["2000::"] #EVERYTHING AFTER THE YEAR 2000
IDX <- index(EXPIRY)
NEXT.EXPIRY <- as.Date("2024-02-22")
IDX <- c(IDX, NEXT.EXPIRY)
# Calculate the mean and stdev for the option prices
MEAN = function(calculateUNTIL) #calculate the mean until a certain date
{
  tmp <- tmp[paste0("::", calculateUNTIL)]
  tmp <- ROC(tmp, type = "discrete")
  tmp[is.na(tmp)]<-0 # if there are any NA's set those to 0
  mean(tmp)
}

STDEV = function(calculateUNTIL) #calculate the mean until a certain date
{
  tmp <- tmp[paste0("::", calculateUNTIL)]
  tmp <- ROC(tmp, type = "discrete")
  tmp[is.na(tmp)]<-0 # if there are any NA's set those to 0
  sd(tmp)
}
means <- do.call(rbind, lapply(as.list(IDX), MEAN))

stdevs <- do.call(rbind, lapply(as.list(IDX), STDEV))

#CALCULATE THE DIFFERENCE BETWEEN THE DATES OF THE OPTIONS EXPIRATION DATES
days = as.numeric(diff(IDX))
#Now we will write the Monte Carlo function
MONTE.CARLO = function(sim, iter, LastIter)
{
  simulations <- sim
  N <- days[iter]
  STK_PRC <- as.numeric(EXPIRY[iter])
  MEAN <- means[iter]
  STDEV <- stdevs[iter]
  stock_prices <- c()
  
  for(i in seq(simulations))
  {
    stock_prices <- c(stock_prices, stk_ret(STK_PRC = STK_PRC, N=N, MEAN=MEAN, STDEV = STDEV))
  }
  
  START <- as.data.frame(round(STK_PRC, 2))
  START.DATE = index(EXPIRY[iter])
  PROBS = as.data.frame(t(round(quantile(stock_prices, probs = seq(0,1,0.05)), 2)))
  
  if(iter == LastIter) #if this is the last iteration
  {
    END <- as.data.frame(NA)
    END.DATE = as.data.frame(NA)
  }else{
    END <- as.data.frame(as.numeric(round(EXPIRY[iter+1], 2)))
    END.DATE = index(EXPIRY[iter+1])
  }
  all <- cbind(START, START.DATE, PROBS, END, END.DATE)
  
  #rename columns
  colnames(all) <- c("START.PRC", "START.DATE", "0%", "5%", "10%", "15%",
                     "20%", "25%", "30%", "35%", "40%", "45%", "50%", "55%",
                     "60%", "65%", "70%", "75%", "80%", "85%", "90%", "95%",
                     "100%", "END.PRC", "END.DATE")
  all

}

p <- pblapply(as.list(1:length(days)), function(x){
  MONTE.CARLO(sim=10000,iter = x, LastIter = length(days))
})

p <- do.call(rbind, p)

plot(p$END.PRC, type = "l")
lines(p$`0%`, col = "red")
lines(p$`100%`, col = "green")

#Number of months
nMo <- nrow(p) #defined by the number of rows

# Number of times it closes above 100% threshold
sum(as.numeric(na.omit(ifelse(p$END.PRC > p$`100%`,1, 0))))/nMo
## [1] 0.08098592
#Number of times it closes below 0%
sum(as.numeric(na.omit(ifelse(p$END.PRC > p$`0%`,1, 0))))/nMo
## [1] 0.9471831