MATLAB Convert R Code to MATLAB for Detrended Cross Correlation Analysis

Click For Summary
The discussion focuses on converting R code for Detrended Cross Correlation Analysis (DCCA) into MATLAB. The original R code calculates beta coefficients and various statistical metrics for two time series. Participants emphasize the importance of understanding the R code and its libraries, such as tseries and fracdiff, to facilitate the conversion process. There is also a suggestion to explore running MATLAB from within R as an alternative approach. Overall, a clear comprehension of the R code structure and functions is deemed essential for successful translation to MATLAB.
Tilfani
Messages
11
Reaction score
0

Homework Statement


I would like to convert this R code into matlab. This code perform a calculation of beta coefficient based on Detrended cross correlation analysis.

Homework Equations


Code: 
require(tseries)

require(fracdiff)

require(matrixStats)DCCA_beta_avg<-function(y,x,smin,smax,step){

  XTX<-var(x)*(length(x)-1)

  betas<-rep(0,(smax-smin)/step+1)

  for(s in seq(smin,smax,by=step)){

    betas[(s-smin)/step+1]<-DCCA_beta_sides(y,x,s)

  }

  DCCA_beta<-mean(betas)

  DCCA_res<-(y-DCCA_beta*x)-mean(y-DCCA_beta*x)

  DCCA_sigma2<-sum(DCCA_res^2)/(length(DCCA_res)-2)

  DCCA_SE<-sqrt(DCCA_sigma2/XTX)

  DCCA_R2<-1-var(DCCA_res)/var(y)

  OLS_beta<-lm(y~x)$coefficients[2]

  OLS_res<-(y-OLS_beta*x)-mean(y-OLS_beta*x)

  OLS_sigma2<-sum(OLS_res^2)/(length(OLS_res)-2)

  OLS_SE<-sqrt(OLS_sigma2/XTX)

  OLS_R2<-1-var(OLS_res)/var(y)

  return(c(OLS_beta,OLS_SE,OLS_R2,DCCA_beta,DCCA_SE,DCCA_R2))

}DCCA_beta<-DCCdccafunction(y,x,s){

  xx<-cumsum(x-mean(x))

  yy<-cumsum(y-mean(y))

  t<-1:length(xx)

  F2sj_xy<-runif(floor(length(xx)/s))

  F2sj_xx<-F2sj_xy

  for(ss in seq(1,(floor(length(xx)/s)*s),by=s)){

    F2sj_xy[(ss-1)/s+1]<-sum((summary(lm(xx[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals)*(summary(lm(yy[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals))/(s-1)

    F2sj_xx[(ss-1)/s+1]<-sum((summary(lm(xx[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals)*(summary(lm(xx[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals))/(s-1)

  }

  beta<-mean(F2sj_xy)/mean(F2sj_xx)

  return(beta)

}DCCA_beta_F<-function(y,x,s){

  xx<-cumsum(x-mean(x))

  yy<-cumsum(y-mean(y))

  t<-1:length(xx)

  F2sj_xy<-runif(floor(length(xx)/s))

  F2sj_xx<-F2sj_xy

  F2sj_yy<-F2sj_xy

  for(ss in seq(1,(floor(length(xx)/s)*s),by=s)){

    F2sj_xy[(ss-1)/s+1]<-sum((summary(lm(xx[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals)*(summary(lm(yy[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals))/(s-1)

    F2sj_xx[(ss-1)/s+1]<-sum((summary(lm(xx[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals)*(summary(lm(xx[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals))/(s-1)

    F2sj_yy[(ss-1)/s+1]<-sum((summary(lm(yy[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals)*(summary(lm(yy[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals))/(s-1)

  }

  beta<-mean(F2sj_xy)/mean(F2sj_xx)

  return(c(beta,mean(F2sj_xx),mean(F2sj_yy)))

  #return(c(beta,sum(F2sj_xx),sum(F2sj_yy)))

}DCCA_beta_SE<-function(y,x,s){

  r<-DCCA_beta_F(y,x,s)

  beta<-r[1]

  yhat<-beta*x

  alpha<-mean(y)-beta*mean(x)

  res<-y-yhat

  residuals<-res-mean(res)

  resres<-cumsum(residuals-mean(residuals))

  F2sj_res<-runif(floor(length(residuals)/s))

  t<-1:length(resres)

  for(ss in seq(1,(floor(length(residuals)/s)*s),by=s)){

    F2sj_res[(ss-1)/s+1]<-sum((summary(lm(resres[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals)*(summary(lm(resres[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals))/(s-1)

  }

  #SE<-mean(residuals^2)/((length(residuals)-2)*r[2])

  SE<-mean(F2sj_res)/((length(residuals)-2)*r[2])

  SE_a<-(mean(F2sj_res)/r[2])*(sum(x^2)/(length(residuals)*(length(residuals)-2)))

  R<-1-mean(F2sj_res)/(r[3])

  return(c(alpha,sqrt(SE_a),beta,sqrt(SE),R))

}DCCA_beta_SE_F<-function(y,x,s){

  r<-DCCA_beta_F(y,x,s)

  beta<-r[1]

  yhat<-beta*x

  alpha<-mean(y)-beta*mean(x)

  res<-y-yhat

  residuals<-res-mean(res)

  res_R<-y-x

  resres<-cumsum(residuals-mean(residuals))

  resres_R<-cumsum(res_R)

  F2sj_res<-runif(floor(length(residuals)/s))

  F2sj_res_R<-runif(floor(length(res_R)/s))

  t<-1:length(resres)

  for(ss in seq(1,(floor(length(residuals)/s)*s),by=s)){

    F2sj_res[(ss-1)/s+1]<-sum((summary(lm(resres[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals)*(summary(lm(resres[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals))/(s-1)

    F2sj_res_R[(ss-1)/s+1]<-sum((summary(lm(resres_R[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals)*(summary(lm(resres_R[ss:(ss+s-1)]~t[ss:(ss+s-1)]))$residuals))/(s-1)

  }

  #SE<-mean(residuals^2)/((length(residuals)-2)*r[2])

  #SE<-mean(F2sj_res)/((length(residuals)-2)*r[2])

  #SE<-mean(F2sj_res)/((length(F2sj_res)-2)*r[2]) #controlling for uncertainty connected to scales (higher scales have higher uncertainty due to lower number of blocks)

  SE<-mean(F2sj_res)/(ceiling(length(residuals)/s)*r[2]) #loosing d.f. due to fitting a and b in each box

  #SE_a<-(mean(F2sj_res)/r[2])*(sum(x^2)/(length(residuals)*(length(residuals)-2)))

  #SE_a<-(mean(F2sj_res)/r[2])*(sum(x^2)/(length(F2sj_res)*(length(F2sj_res)-2))) #controlling for uncertainty connected to scales (higher scales have higher uncertainty due to lower number of blocks)

  SE_a<-(mean(F2sj_res)/r[2])*(sum(x^2)/(length(residuals)*ceiling(length(residuals)/s))) #loosing d.f. due to fitting a and b in each box

  R<-1-mean(F2sj_res)/(r[3])

  #SSR_U<-sum(residuals^2)

  SSR_U<-sum(F2sj_res)

  #SSR_R<-sum((y-x)^2) #specific null: alpha=0, beta=1

  SSR_R<-sum(F2sj_res_R)

  #F_stat<-((SSR_R-SSR_U)/(SSR_U))*((length(residuals)-2)/2)

  #F_stat<-((SSR_R-SSR_U)/(SSR_U))*((length(F2sj_res)-2)/2) #controlling for uncertainty connected to scales (higher scales have higher uncertainty due to lower number of blocks)

  F_stat<-((SSR_R-SSR_U)/(SSR_U))*(ceiling(length(residuals)/s)/2) #loosing d.f. due to fitting a and b in each box

  F_p<-pf(F_stat,2,length(F2sj_res)-2,lower.tail=FALSE)

  return(c(alpha,sqrt(SE_a),beta,sqrt(SE),R,F_stat,F_p))

}DCCA_beta_s<-function(y,x,smin,smax,step){

  results<-matrix(rep(0,6*((smax-smin)/step+1)),ncol=6)

  for(s in seq(smin,smax,by=step)){

    beta<-DCCA_beta_SE(y,x,s)

    results[((s-smin)/step+1),1]<-s

    results[((s-smin)/step+1),2]<-beta[1]

    results[((s-smin)/step+1),3]<-beta[2]

    results[((s-smin)/step+1),4]<-beta[3]

    results[((s-smin)/step+1),5]<-beta[4]

    results[((s-smin)/step+1),6]<-beta[5]

  }

  return(results)

}DCCA_beta_s_F<-function(y,x,smin,smax,step){

  results<-matrix(rep(0,10*((smax-smin)/step+2)),ncol=10)

  for(s in seq(smin,smax,by=step)){

    beta<-DCCA_beta_SE_F(y,x,s)

    results[((s-smin)/step+1),1]<-s

    results[((s-smin)/step+1),2]<-beta[1]

    results[((s-smin)/step+1),3]<-beta[2]

    results[((s-smin)/step+1),4]<-2*pnorm(abs(beta[1]/beta[2]),lower.tail=FALSE)#p-value for null=0

    results[((s-smin)/step+1),5]<-beta[3]

    results[((s-smin)/step+1),6]<-beta[4]

    results[((s-smin)/step+1),7]<-2*pnorm(abs((beta[3]-1)/beta[4]),lower.tail=FALSE)#p-value for null=1

    results[((s-smin)/step+1),8]<-beta[5]

    results[((s-smin)/step+1),9]<-beta[6]

    results[((s-smin)/step+1),10]<-beta[7]

  }

  #results[(smax-smin)/step+2,2]<-mean(results[1:(dim(results)[1]-1),2])#A

  #results[(smax-smin)/step+2,5]<-mean(results[1:(dim(results)[1]-1),5])#B

  results[(smax-smin)/step+2,2]<-sum(results[1:(dim(results)[1]-1),2]*results[1:(dim(results)[1]-1),8])/sum(results[1:(dim(results)[1]-1),8])#A as R2(s) weighted

  results[(smax-smin)/step+2,5]<-sum(results[1:(dim(results)[1]-1),5]*results[1:(dim(results)[1]-1),8])/sum(results[1:(dim(results)[1]-1),8])#B as R2(s) weighted

  results[(smax-smin)/step+2,3]<-sqrt((sum(x^2)/length(x))*sum((y-results[(smax-smin)/step+2,2]-results[(smax-smin)/step+2,5]*x)^2)/((length(y)-dim(results)[1]+1)*sum((x-mean(x))^2)))#SE_A

  results[(smax-smin)/step+2,4]<-2*pnorm(abs(results[(smax-smin)/step+2,2]/results[(smax-smin)/step+2,3]),lower.tail=FALSE)#p-value for null=0

  results[(smax-smin)/step+2,6]<-sqrt(sum((y-results[(smax-smin)/step+2,2]-results[(smax-smin)/step+2,5]*x)^2)/((length(y)-dim(results)[1]+1)*sum((x-mean(x))^2)))#SE_B

  results[(smax-smin)/step+2,7]<-2*pnorm(abs((results[(smax-smin)/step+2,5]-1)/results[(smax-smin)/step+2,6]),lower.tail=FALSE)#p-value for null=1

  results[(smax-smin)/step+2,8]<-1-sum((y-results[(smax-smin)/step+2,2]-results[(smax-smin)/step+2,5]*x)^2)/sum(y^2)#R2

  results[(smax-smin)/step+2,9]<-((length(x)-2)/2)*((results[(smax-smin)/step+2,8]-(1-(sum((y-x)^2))/(sum((y-mean(y))^2))))/(1-results[(smax-smin)/step+2,8]))#F_test_R2_based

  results[(smax-smin)/step+2,10]<-pf(results[(smax-smin)/step+2,9],2,length(y)-2,lower.tail=FALSE)#F_test p_val

  return(results)

}

The Attempt at a Solution



I would like to obtain a MATLAB code.[/B]
 
Last edited by a moderator:
Physics news on Phys.org
I added code tags for readability -- please always use code tags to post programs here. Thanks.

So this is for schoolwork? Can you show us your Attempt at the Solution? We don't do your schoolwork for you here.
 
berkeman said:
I added code tags for readability -- please always use code tags to post programs here. Thanks.

So this is for schoolwork? Can you show us your Attempt at the Solution? We don't do your schoolwork for you here.
It's not a schoolwork.
I would like to perform an anlysis of two time series x and y, my aim is the perform a particular regression approach, on them.
 
Okay, let me move it to the Programming forum then...

ADD -- Have you started the MATLAB code? Can you post what you have done so far?
 
berkeman said:
Okay, let me move it to the Programming forum then...

ADD -- Have you started the MATLAB code? Can you post what you have done so far?
Not yet, i didn't understant what this codes do, i don't understand R programming language.
 
Tilfani said:
Not yet, i didn't understant what this codes do, i don't understand R programming language.
If you don't understand the R code, you will not be able to convert it to MATLAB code. If I were you I would start by looking at the R documentation.

The three lines below seem to be libraries that the code uses. If so, it would be useful to look at these libraries to see what the functions in them do, what parameters are required, and what kinds of values they return.
Code: 
require(tseries)
require(fracdiff)
require(matrixStats)

It would also be useful to understand what this notation means. I don't know R, so I can't offer any advice, but I would be looking to see what DCCA_beta_avg is, and what the notation <- means.
Code: 
DCCA_beta_avg<-function(y,x,smin,smax,step)
From the line below, it's possible that R uses <- for assignment operations, but I can't say for sure.
Code: 
XTX<-var(x)*(length(x)-1)
 
Step 1 is to find the person who wrote the R code without a single comment and poke him in the eye. Do you even know what the input and the output looks like? If so, you may be able to write MATLAB code without knowing every detail of the R code.
 
  • Like
Likes berkeman
Why do you want to convert it to matlab? It is possible to run MATLAB from inside R and use the results for computation. This would be probably easier.
 
  • Like
Likes berkeman
Mark44 said:
Code: 
require(tseries)
require(fracdiff)
require(matrixStats)
Yes, they can be found on the CRAN archive, including their descriptions.
 

Similar threads

MATLAB Problem with Runge-Kutta script file
  • · Replies 4 ·
Replies
4
Views
1K
MATLAB Transforming part of matlab code to Fortran90
  • · Replies 8 ·
Replies
8
Views
2K
MATLAB How do I write this chua oscillator to an audio file in matlab?
  • · Replies 5 ·
Replies
5
Views
3K
MATLAB How to plot flow data from an orifice in MATLAB without overlapping the plots?
  • · Replies 2 ·
Replies
2
Views
2K
MATLAB To what extent is the fit to experimental data good?
  • · Replies 12 ·
Replies
12
Views
4K
MATLAB Fixing a Matlab Code: Help Needed for Fourier Transform Calculation"
  • · Replies 29 ·
Replies
29
Views
5K
MATLAB What is the issue with my step function code and how can I fix it?
  • · Replies 3 ·
Replies
3
Views
2K
MATLAB Crank-Nicholson in 2D with MATLAB
  • · Replies 1 ·
Replies
1
Views
4K
MATLAB System of differential equations
  • · Replies 9 ·
Replies
9
Views
2K
Mathematica How Can I Solve Complex PDEs Using Mathematica's NDSolve?
  • · Replies 2 ·
Replies
2
Views
3K