How to obtain generalizations in R package Biwavelet?

Question

I’m preparing a coursework devoted to diversification strategies for stock market. To measure interdependence of national stock markets I use R with Biwavelet package.

While graphs are really good for noticing connections in various periods, scientific work requires numerical results. How can I obtain numbers that will in short generalise the graphs?

The data is here:

https://docs.google.com/spreadsheets/d/1i5Yvufe2IUrx42sW9WR0hKMeCTyNyrSGn-fpenE2obQ/edit?usp=sharing

library(biwavelet)

wtc_Modfd <- function (d1, d2, pad = TRUE, dj = 1/12, s0 = 2 * dt, J1 = NULL, 
                       max.scale = NULL, mother = "morlet", param = -1, lag1 = NULL, 
                       sig.level = 0.95, sig.test = 0, nrands = 300, quiet = FALSE) 
{
  mother <- match.arg(tolower(mother), MOTHERS)
  checked <- check.data(y = d1, x1 = d2)
  xaxis <- d1[, 1]
  dt <- checked$y$dt
  t <- checked$y$t
  n <- checked$y$n.obs
  if (is.null(J1)) {
    if (is.null(max.scale)) {
      max.scale <- (n * 0.17) * 2 * dt
    }
    J1 <- round(log2(max.scale/s0)/dj)
  }
  if (is.null(lag1)) {
    d1.ar1 <- arima(d1[, 2], order = c(1, 0, 0),method="ML")$coef[1] # Added method to be ML to avoid error: non-stationary AR part from CSS
    d2.ar1 <- arima(d2[, 2], order = c(1, 0, 0),method="ML")$coef[1] # Added method to be ML to avoid error: non-stationary AR part from CSS
    lag1 <- c(d1.ar1, d2.ar1)
  }
  wt1 <- wt(d = d1, pad = pad, dj = dj, s0 = s0, J1 = J1, max.scale = max.scale, 
            mother = mother, param = param, sig.level = sig.level, 
            sig.test = sig.test, lag1 = lag1[1])
  wt2 <- wt(d = d2, pad = pad, dj = dj, s0 = s0, J1 = J1, max.scale = max.scale, 
            mother = mother, param = param, sig.level = sig.level, 
            sig.test = sig.test, lag1 = lag1[2])
  d1.sigma <- sd(d1[, 2], na.rm = T)
  d2.sigma <- sd(d2[, 2], na.rm = T)
  s.inv <- 1/t(wt1$scale)
  s.inv <- matrix(rep(s.inv, n), nrow = NROW(wt1$wave))
  smooth.wt1 <- smooth.wavelet(s.inv * (abs(wt1$wave)^2), dt, 
                               dj, wt1$scale)
  smooth.wt2 <- smooth.wavelet(s.inv * (abs(wt2$wave)^2), dt, 
                               dj, wt2$scale)
  coi <- pmin(wt1$coi, wt2$coi, na.rm = T)
  CW <- wt1$wave * Conj(wt2$wave)
  CW.corr <- (wt1$wave * Conj(wt2$wave) * max(wt1$period))/matrix(rep(wt1$period, 
                                                                      length(t)), nrow = NROW(wt1$period))
  power <- abs(CW)^2
  power.corr <- (abs(CW)^2 * max.scale)/matrix(rep(wt1$period, 
                                                   length(t)), nrow = NROW(wt1$period))
  smooth.CW <- smooth.wavelet(s.inv * (CW), dt, dj, wt1$scale)
  rsq <- abs(smooth.CW)^2/(smooth.wt1 * smooth.wt2)
  phase <- atan2(Im(CW), Re(CW))
  if (nrands > 0) {
    signif <- wtc.sig(nrands = nrands, lag1 = lag1, dt = dt, 
                      ntimesteps = n, pad = pad, dj = dj, J1 = J1, s0 = s0, 
                      max.scale = max.scale, mother = mother, sig.level = sig.level, 
                      quiet = quiet)
  }
  else {
    signif <- NA
  }
  results <- list(coi = coi, wave = CW, wave.corr = CW.corr, 
                  power = power, power.corr = power.corr, rsq = rsq, phase = phase, 
                  period = wt1$period, scale = wt1$scale, dt = dt, t = t, 
                  xaxis = xaxis, s0 = s0, dj = dj, d1.sigma = d1.sigma, 
                  d2.sigma = d2.sigma, mother = mother, type = "wtc", 
                  signif = signif)
  class(results) <- "biwavelet"
  return(results)
}

US <- read.csv('US.txt', header = FALSE, sep = ",", dec = ".")[c('V3','V5')]
US <- US[which(US['V3'] == '20071025'):nrow(US),]

India <- read.csv('India.txt', header = FALSE, sep = ",", dec = ".")[c('V3','V5')]
India <- India[which(India['V3'] == '20071025'):nrow(India),]

joint <- merge(India, US, by=c('V3'))
joint <- na.omit(joint)

nrands=10 #Number of iterations. Higher is better (>500)

# Define two sets of variables with time stamps
DATE=1:nrow(joint)

WA = cbind(DATE, joint['V5.x'])
WB = cbind(DATE, joint['V5.y'])

# Specify the number of iterations. The more, the better (>1000).  For the
# purpose of this tutorial, we just set it = 10
nrands = 100

#Main calculations
wtc.AB = wtc(WA, WB, nrands = nrands)

# wt1 = wtc_Modfd(WA, WB, nrands = nrands)

# Plotting a graph
par(oma = c(0, 0, 0, 1), mar = c(5, 4, 5, 5) + 0.1)
plot(wtc.AB, plot.phase = T, lty.coi = 1, col.coi = "grey", lwd.coi = 2, 
     lwd.sig = 1, arrow.lwd = 0.03, arrow.len = 0.12, ylab = "Period (days)", xlab = "Time (days)", 
     plot.cb = TRUE, main = "")

# Adding grid lines
n = length(WA[, 1])
abline(v = seq(250, n, 250), h = 1:4, col = "brown", lty = 1, lwd = 1)
# Defining x labels
axis(side = 3, at = c(44,320,596,838,1078,1324,1572,1820,2060,2302,2541,2780,3022,3266,3508), 
     labels = c('2008','2009', '2010','2011','2012','2013','2014','2015',
                                                  '2016','2017','2018','2019','2020','2021','2022'))

As for the result i want a table with average covariation for various periods (4,12,16,32,64,128,256 days) measured for the whole time and particular time spans.