conert R code into matlab : I would like to convert this code from R into matlab.
this aimed to perform a calculations of regression parameters a detrended cross correlation analysis.
there is a way to do it automatically, or we should do it line by line.
This analysis should be performed on tiem series.
x and y denotes two time series with same length.

   

     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)
    }


