我对绝对利率假设的数值检验

 # clear environment rm(list = ls()); gc() ## load libs library(data.table) library(ggplot2) library(magrittr) ## set WD -------------------------------- # your dir here ... ## set vars ------------- currs <- c( 'usd', 'eur', 'chf', 'gbp' ) ############ ## RUN SIMULATION LOOP ------------------------------- simuls <- 1000L simul_dt <- data.table() for( s in seq_len(simuls) ) { #set.seed(111) usd <- runif(1) eur <- runif(1) chf <- runif(1) gbp <- runif(1) # snapshot of values at time t eurusd <- 1.12012 gbpusd <- 1.30890 eurchf <- 1.14135 eurgbp <- 0.85570 gbpchf <- 1.33373 usdchf <- 1.01896 ## symbolic task ------------ express <- expression( (eurusd - eur / usd) ^ 2 + (gbpusd - gbp / usd) ^ 2 + (eurchf - eur / chf) ^ 2 + (eurgbp - eur / gbp) ^ 2 + (gbpchf - gbp / chf) ^ 2 + (usdchf - usd / chf) ^ 2 ) ## define gradient and iterate to make descent to zero -------------- iter_max <- 1e+3 lr <- 1e-3 min_tolerance <- 0.00001 rm(grad_desc_func) grad_desc_func <- function( lr, curr_list ) { derivs <- character(length(curr_list)) deriv_vals <- numeric(length(curr_list)) grads <- numeric(length(curr_list)) # symbolic derivatives derivs <- sapply( curr_list, function(x){ D(express, x) } ) # derivative values deriv_vals <- sapply( derivs, function(x){ eval(x) } ) # gradient change values -deriv_vals * lr } ## get gradient values ---------- progress_list <- list() for( i in seq_len(iter_max) ) { grad_deltas <- grad_desc_func(lr, curr_list = currs) currency_vals <- sapply( currs , function(x) { # update currency values current_val <- get(x, envir = .GlobalEnv) new_delta <- grad_deltas[x] if(new_delta > -1 & new_delta < 1) { new_delta = new_delta } else { new_delta = sign(new_delta) } new_val <- current_val + new_delta if(new_val > 0 & new_val < 2) { new_val = new_val } else { new_val = current_val } names(new_val) <- NULL # change values of currencies by gradient descent step in global env assign(x, new_val , envir = .GlobalEnv) # save history of values for later plotting new_val } ) progress_list[[i]] <- c( currency_vals, eval(express) ) if( eval(express) < min_tolerance ) { break('solution was found') } } ## check results ---------- # print( # paste0( # 'Final error: ' # , round(eval(express), 5) # ) # ) # # print( # round(unlist(mget(currs)), 5) # ) progress_dt <- rbindlist( lapply( progress_list , function(x) { as.data.frame(t(x)) } ) ) colnames(progress_dt)[length(colnames(progress_dt))] <- 'error' progress_dt[, steps := 1:nrow(progress_dt)] progress_dt_melt <- melt( progress_dt , id.vars = 'steps' , measure.vars = colnames(progress_dt)[colnames(progress_dt) != 'steps'] ) progress_dt_melt[, simul := s] simul_dt <- rbind( simul_dt , progress_dt_melt ) } ggplot(data = simul_dt) + facet_wrap(~ variable, scales = 'free') + geom_line( aes( x = steps , y = value , group = simul , color = simul ) ) + scale_y_log10() + theme_minimal()